Jordan形$O(n^4)$算法
Computing Jordan canonical form
Definition 1 root subspace for \(\lambda\) is the subspace satisfying \(\exist k, (\mathcal{A}-\lambda\mathcal{I})^k\alpha=0\)
Lemma 1 The sum of different root subspace is direst sum
proof
\(f(x)\) is the characteristic polynomial's factor of \(\lambda_1\) , \(g(x)\) is the remaining part
Now \(f(x)\) is annihilator for \(W_{\lambda _1}\) because the whole characteristic polynomial is annihilator, but \(g(x)\) cannot annihilate \(W_{\lambda _1}\) (consider the diagonal of the matrices)
we have \((f(x),g(x))=1\) now we claim that \(g(\mathcal{A})\) is invertible. There exist \(u(x)g(x)=1 \bmod f(x)\)
so \(u(\mathcal{A})g(\mathcal{A})=\mathcal{I}\)
to prove \(W_{\lambda_1}\oplus W_{\lambda_2}\dots\oplus W_{\lambda_t}\) we only need to prove vector \(0\) have only one representation.
say \(0=u_1+u_2+\dots+u_t\) . applying \(g(\mathcal{A})\) to both side yield \(g(\mathcal{A})u_1=0\) . Then applying \(u(\mathcal{A})\) gives the result. \(\square\)
Definition 2 cyclic subspace \(C(\alpha)\) is the subspace generated by \(\alpha \ ,\mathcal{A}\alpha,\ \mathcal{A}^2\alpha \dots\) and denote the last none zero term of the sequence \(last(\alpha)\)
(it's easy to prove \(\alpha \ ,\mathcal{A}\alpha,\ \mathcal{A}^2\alpha \dots\) are independent.
Lemma 2 Nilpotent operator is the direct sum of some cyclic subspace
proof
we give a constructive and algorithmic proof.
Let the degree of minimal polynomial be \(d\).
First start with solving \(AX=0\) ,then \(A^2X=0\) ..... \(A^dX=0\) The solution space is expanding. Denote \(X_i\) the set of basis that is newly inserted after solving the i-th equation.
Denote the vectors in \(X_d\) \(\alpha_1,\alpha_2\dots\alpha_k\) , in \(X_{d-1}\) \(\beta_1,\beta_2\dots\beta_l\) in \(X_{d-2}\) \(\gamma_1,\gamma_2\dots\gamma_m\)
First we claim that \(C(\alpha_1)\oplus C(\alpha_2)\dots \oplus C(\alpha_k)\) .
applying \(\mathcal{A}^{d-1}\) to both side yields
now
is in \(Ker(\mathcal{A}^{d-1})\) ,but \(\alpha_i\) are newly inserted, if this is not \(0\), some \(\alpha\) can be represented by the basis of \(Ker(\mathcal{A}^{d-1})\) ,a contradiction.
so \(\forall i,\mu_{i0}=0\) because \(\alpha_i\) are linearly independent.
now apply \(\mathcal{A}^{d-2}\) and similarly \(\forall i,\mu_{i1}=0\). the remaining proof is obvious.
At that time, we want to insert \(C(\beta_i)\) but not all of them should be inserted.
Firstly I thought if \(\beta\) is linearly independent with all \(\mathcal{A}\alpha\) , \(C(\beta)\) and \(C(\alpha)\) is direct sum. but that was wrong.
consider the following example \(A=\)
0 0 1 2 3
0 0 0 1 -1
0 0 0 1 2
0 0 0 0 0
0 0 0 0 0
\(X_1\) is
1 0
0 1
0 0
0 0
0 0
\(X_2\) is
0 0
0 0
1 0
0 -2
0 1
\(X_3\) is
0
0
0
0
1
now \(\alpha\) \(\mathcal{A} \alpha\) \(\mathcal{A}^2 \alpha\) are
0 3 2
0 -1 0
0 2 0
0 0 0
1 0 0
now choose in \(X_2\) \(\beta =\)
0
0
1
0
0
which is linearly independent with \(\mathcal{A}\alpha\) but \(\mathcal{A}\beta=\)
1
0
0
0
0
is linearly independent with \(\mathcal{A}^2 \alpha\)
so it occurred to me that I should maintain a linearly independent kernel.
apply \(\mathcal{A}^{d-1}\) we have.
now apply \(\mathcal{A}^{d-2}\)
because the \(\beta\) we choose (\(\beta\) here is not all the \(\beta\) )to add here satisfy \(last(\beta)\) is linearly independent with \(last(\alpha)\)
so \(\forall i, \mu_{i,1}=\lambda_{i,0}=0\) continue the process, we can prove we get some direct sum of cyclic subspace.
but is the sum the complete space \(V\) ?
Note that originally \(X_1, X_2 \dots\) are a basis.
when we come to the k-th layer(corresponding to inserting \(X_k\)) , we can consider this algorithm as substituting some vectors in \(X_k\) with the previously inserted vectors. For example, use \(\mathcal{A} \alpha\) to change out some \(\beta\) . Because originally \(last(X_{d-1})\) is linearly independent with size \(|X_{d-1}|=|last(X_{d-1})|\) , finally the vectors in the answer set form a set with size \(|last(X_{d-1})|\) too.(the basis of a space is a fixed number).
So finally we will get a set with size \(|X_{d-1}|\) in place of \(X_{d-1}\). For \(X_k\) it is similar. \(\square\)
Theorem There is a basis under which the matrix is in Jordan form.
proof
For every eigenvalue, run algorithm of lemma 2. And the union of vectors is a Jordan basis.
Furthermore, it is easy to see that the number of different size of Jordan square is determined by the solution of \((A-\lambda^kI)X=0\) so it is determined by \(rank(A-\lambda I)^k\) . Which is invariant for similar relation.
So Jordan form is canonical form of similar matrices. \(\square\)
Corollary Jordan form and its transfer matrix can be computed in \(O(n^4)\).
proof
Computing the solution for \(A^kX=0\) need \(Cn^4\) operations if implemented using ordinary method.
Computing the chain of a vector \(\alpha\) (the space \(C(\alpha)\)) need \(C\sum d_i\times n^2=Cn^3\). Maintaining the kernel need \(C\sum d_i^3=O(n^3)\)
So the complexity is \(O(n^4)\) \(\square\)
This method is a clear and straightforward geometry method. It is easy to implement compared to algebraic method of \(\lambda\)-matrix. And the proof of lemma 2 is different from classical textbook because I want to prove the correctness of the algorithm.
There are better algorithms for computing Jordan form. For example, An \(O(n^3)\) Algorithm for the Frobenius Normal Form Storjohann 1998 is an easy faster way.
I did not find a fast algorithm that can compute transfer matrix of Jordan form online.
but use idea of some reduction. It may be done with better complexity.
The code below needs eigenvalues.
#include<bits/stdc++.h>
#define de fprintf(stderr,"on-%d\n",__LINE__)
#define pb push_back
#define fr(i,a,b) for(int i(a),end##i(b);i<=end##i;i++)
#define fd(i,a,b) for(int i(a),end##i(b);i>=end##i;i--)
#define reduce Reduce
#define apply Apply
using namespace std;
const int N=30;
int m;
struct frac{
int q,p;
frac(){
q=0;p=1;
}
void out(){
if(p==1)printf("%7d ",q);
else printf("%3d/%3d ",q,p);
}
};
frac zero;
frac make(int x,int y){
frac A;
A.q=x;A.p=y;
return A;
}
frac r[N];
int d[N];
frac V(int x){
return make(x,1);
}
frac read(){
int x;scanf("%d",&x);
return V(x);
}
int gcd(int x,int y){
if(!x)return y;
return gcd(y%x,x);
}
frac reduce(frac &A){
int g=gcd(A.p,A.q);
A.p/=g;A.q/=g;
if(A.p<0)A.p=-A.p,A.q=-A.q;
return A;
}
frac reduce(frac A){
int g=gcd(A.p,A.q);
A.p/=g;A.q/=g;
if(A.p<0)A.p=-A.p,A.q=-A.q;
return A;
}
frac operator +(const frac &A,const frac &B){
return reduce(make(A.q*B.p+B.q*A.p,A.p*B.p));
}
frac operator -(const frac &A,const frac &B){
return reduce(make(A.q*B.p-B.q*A.p,A.p*B.p));
}
frac operator *(const frac &A,const frac &B){
assert(A.p);
return reduce(make(A.q*B.q,A.p*B.p));
}
frac operator /(const frac &A,const frac &B){
return reduce(make(A.q*B.p,A.p*B.q));
}
struct matrix{
frac A[N][N];
frac* operator [](int x){return A[x];}
void out(){
fr(i,1,m){
fr(j,1,m)A[i][j].out();
printf("\n");
}
}
};
matrix I;
void init(){
fr(i,1,m)fr(j,1,m)I[i][j]=V(i==j);
}
matrix operator *(matrix A,matrix B){
matrix C;
fr(i,1,m)fr(j,1,m){
C[i][j]=V(0);
fr(k,1,m)C[i][j]=C[i][j]+A[i][k]*B[k][j];
}
return C;
}
matrix operator *(frac v,matrix A){
fr(i,1,m)fr(j,1,m)A[i][j]=v*A[i][j];
return A;
}
matrix operator +(matrix A,matrix B){
fr(i,1,m)fr(j,1,m)A[i][j]=A[i][j]+B[i][j];
return A;
}
matrix operator -(matrix A,matrix B){
fr(i,1,m)fr(j,1,m)A[i][j]=A[i][j]-B[i][j];
return A;
}
vector<frac> apply(matrix &A,const vector<frac>& B){
vector<frac> C;C.resize(m+1);
fr(i,1,m)fr(k,1,m){
C[i]=C[i]+A[i][k]*B[k];
}
return C;
}
bool p[N];
int to[N];
void rowreduce(matrix &A,int n,int m,int op=0){
fr(i,1,m)to[i]=p[i]=0;
fr(i,1,n){
int s=0;
fr(j,1,m)if(!p[j]&&A[i][j].q){
s=j;break;
}
if(!s)continue;
to[i]=s;
p[s]=1;
fr(j,op*i+1,n)if(j!=i&&A[j][s].q){
frac t=A[j][s]/A[i][s];
fr(k,1,m)A[j][k]=A[j][k]-t*A[i][k];
}
}
}
void swap(frac* A,frac *B){
fr(i,1,m)swap(A[i],B[i]);
}
matrix inverse(matrix A){
matrix B=I;
fr(i,1,m){
int s=0;
fr(j,i,m)if(A[j][i].q){
s=j;break;
}
assert(s);
swap(A[i],A[s]);swap(B[i],B[s]);
const frac t=A[i][i];
fr(j,1,m){
A[i][j]=A[i][j]/t;
B[i][j]=B[i][j]/t;
}
fr(j,1,m)if(j!=i&&A[j][i].q){
const frac t=A[j][i];
fr(k,1,m){
A[j][k]=A[j][k]-t*A[i][k];
B[j][k]=B[j][k]-t*B[i][k];
}
}
}
return B;
}
vector<vector<frac> > solve(matrix A){
vector<vector<frac> >W;
rowreduce(A,m,m,0);
fr(i,1,m)if(!p[i]){
vector<frac> T;T.resize(m+1);
T[i]=V(1);
fr(j,1,m)if(to[j]&&A[j][i].q){
T[to[j]]=zero-A[j][i]/A[j][to[j]];
}
W.push_back(T);
}
return W;
}
bool noempty(frac* A){
fr(i,1,m)if(A[i].q)return 1;
return 0;
}
int main(){
#ifdef pig
freopen("pig.in","r",stdin);
freopen("pig.out","w",stdout);
#endif
cin>>m;
init();
matrix A;
fr(i,1,m)fr(j,1,m)A[i][j]=read();
int s;cin>>s;
fr(i,1,s){
r[i]=read();
scanf("%d",&d[i]);
}
matrix Q;int top=0;
fr(i,1,s){
matrix B=A-r[i]*I,C=B;
vector<vector<vector<frac> > > V;
vector<vector<frac> > P;
fr(j,1,d[i]){
vector<vector<frac> > W=solve(C);
matrix A;
fr(k,0,P.size()-1){
fr(l,1,m)A[k+1][l]=P[k][l];
}
fr(k,0,W.size()-1){
fr(l,1,m){
A[k+P.size()+1][l]=W[k][l];
}
}
rowreduce(A,P.size()+W.size(),m,1);
vector<vector<frac> >T;
fr(k,P.size()+1,W.size()+P.size())if(noempty(A[k])){
vector<frac> TT;TT.resize(m+1);
fr(l,1,m)TT[l]=A[k][l];
T.pb(TT);
}
V.pb(T);
P=W;
C=C*B;
}
vector<vector<frac> >ans;
matrix now;int tot=0;
fd(j,V.size()-1,0){
fr(k,0,V[j].size()-1){
vector<frac> T; T.resize(m+1);
fr(l,1,m){
T[l]=V[j][k][l];
}
fr(l,0,j-1)T=apply(B,T);
tot++;
fr(l,1,m)now[tot][l]=T[l];
}
rowreduce(now,tot,m,1);
fr(k,tot-V[j].size()+1,tot)if(noempty(now[k])){
vector<frac> T;T.resize(m+1);
fr(l,1,m)T[l]=V[j][k-(tot-V[j].size()+1)][l];
fr(t,0,j){
ans.pb(T);
T=apply(B,T);
}
}
}
fr(j,0,ans.size()-1){
top++;
fr(k,1,m)Q[k][top]=ans[j][k];
}
}
assert(top==m);
Q.out();
printf("\n");
(inverse(Q)*A*Q).out();
return 0;
}