矩阵的QR分解

  1 #include <cstdio>
  2 #include <cstdlib>
  3 #include <algorithm>
  4 #include <cmath>
  5 #include <cassert>
  6 #include <ctime>
  7 
  8 class MclVector
  9 {
 10 public:
 11     int n;
 12     double *Mat;
 13     /**
 14       type=0: 列向量
 15       type=1: 行向量
 16     **/
 17     int type;
 18 
 19     MclVector() { Mat=NULL; n=0; }
 20     MclVector(int len,double initVal=0.0)
 21     {
 22         n=len;
 23         Mat=new double[n+1];
 24         for(int i=0;i<=n;i++) Mat[i]=initVal;
 25         type=0;
 26     }
 27 
 28 
 29     double operator[](int id) const
 30     {
 31         return Mat[id];
 32     }
 33 
 34 
 35     double& operator[](int id)
 36     {
 37         return Mat[id];
 38     }
 39 
 40     double length() const
 41     {
 42         double sum=0;
 43         for(int i=1;i<=n;i++) sum+=Mat[i]*Mat[i];
 44         return sqrt(sum);
 45     }
 46 
 47     MclVector operator*(double val) const
 48     {
 49         MclVector ans=MclVector(n);
 50         for(int i=1;i<=n;i++) ans[i]=Mat[i]*val;
 51         return ans;
 52     }
 53 
 54     MclVector operator/(double val) const
 55     {
 56         MclVector ans=MclVector(n);
 57         for(int i=1;i<=n;i++) ans[i]=Mat[i]/val;
 58         return ans;
 59     }
 60 
 61     MclVector operator+(const MclVector &newVector) const
 62     {
 63         MclVector ans=MclVector(n);
 64         for(int i=1;i<=n;i++) ans[i]=Mat[i]+newVector[i];
 65         return ans;
 66     }
 67 
 68 
 69     MclVector operator-(const MclVector &newVector) const
 70     {
 71         MclVector ans=MclVector(n);
 72         for(int i=1;i<=n;i++) ans[i]=Mat[i]-newVector[i];
 73         return ans;
 74     }
 75 
 76     MclVector operator*=(double val)
 77     {
 78         for(int i=1;i<=n;i++) Mat[i]=Mat[i]*val;
 79         return *this;
 80     }
 81 
 82     MclVector operator/=(double val)
 83     {
 84         for(int i=1;i<=n;i++) Mat[i]=Mat[i]/val;
 85         return *this;
 86     }
 87 
 88 
 89     MclVector operator+=(const MclVector &newVector)
 90     {
 91         for(int i=1;i<=n;i++) Mat[i]+=newVector[i];
 92         return *this;
 93     }
 94 
 95     MclVector operator-=(const MclVector &newVector)
 96     {
 97         for(int i=1;i<=n;i++) Mat[i]-=newVector[i];
 98         return *this;
 99     }
100 
101     MclVector GetTranspose() const
102     {
103         MclVector ans=*this;
104         ans.type=1;
105         return ans;
106     }
107 
108 
109     void print() const
110     {
111         for(int i=1;i<=n;i++) printf("%8.3lf ",Mat[i]);
112         puts("");
113     }
114 
115 
116 };
117 
118 class MclMatrix
119 {
120 public:
121     int row,col;
122     MclVector *Mat;
123 
124     MclMatrix() {Mat=NULL;}
125     MclMatrix(int _row,int _col,double initVal=0.0)
126     {
127         row=_row;
128         col=_col;
129         Mat=new MclVector[row+1];
130         for(int i=0;i<=row;i++) Mat[i]=MclVector(col,initVal);
131     }
132 
133     void setIdentityMatrix()
134     {
135         for(int i=1;i<=row;i++)
136         {
137             for(int j=1;j<=col;j++)
138             {
139                 if(i==j) Mat[i][j]=1;
140                 else Mat[i][j]=0;
141             }
142         }
143     }
144 
145     MclMatrix GetTranspose() const
146     {
147         MclMatrix ans=MclMatrix(col,row);
148         for(int i=1;i<=ans.row;i++)
149         {
150             for(int j=1;j<=ans.col;j++)
151             {
152                 ans[i][j]=Mat[j][i];
153             }
154         }
155         return ans;
156     }
157 
158     void print() const
159     {
160         for(int i=1;i<=row;i++) Mat[i].print();
161         puts("");
162     }
163 
164     MclVector& operator[](int id) const
165     {
166         return Mat[id];
167     }
168 
169 
170     MclVector& operator[](int id)
171     {
172         return Mat[id];
173     }
174 
175     MclMatrix operator*(const MclMatrix &Matrix) const
176     {
177         MclMatrix ans=MclMatrix(row,Matrix.col);
178         for(int i=1;i<=row;i++)
179         {
180             for(int j=1;j<=Matrix.col;j++)
181             {
182                 for(int k=1;k<=col;k++)
183                 {
184                     ans[i][j]+=Mat[i][k]*Matrix[k][j];
185                 }
186             }
187         }
188         return ans;
189     }
190 
191     MclMatrix operator+(const MclMatrix &Matrix) const
192     {
193         MclMatrix ans=MclMatrix(row,Matrix.col);
194         for(int i=1;i<=row;i++)
195         {
196             for(int j=1;j<=Matrix.col;j++)
197             {
198                 ans[i][j]=Mat[i][j]+Matrix[i][j];
199             }
200         }
201         return ans;
202     }
203 
204     MclMatrix operator-(const MclMatrix &Matrix) const
205     {
206         MclMatrix ans=MclMatrix(row,Matrix.col);
207         for(int i=1;i<=row;i++)
208         {
209             for(int j=1;j<=Matrix.col;j++)
210             {
211                 ans[i][j]=Mat[i][j]-Matrix[i][j];
212             }
213         }
214         return ans;
215     }
216 
217     MclVector GetCol(int colId) const
218     {
219         MclVector ans=MclVector(row);
220         for(int i=1;i<=row;i++) ans[i]=Mat[i][colId];
221         return ans;
222     }
223     MclVector GetRow(int rowId) const
224     {
225         MclVector ans=MclVector(row);
226         for(int i=1;i<=col;i++) ans[i]=Mat[rowId][i];
227         return ans;
228     }
229 
230 
231     MclMatrix operator*=(const MclMatrix &Matrix)
232     {
233         return *this=*this*Matrix;
234     }
235     MclMatrix operator+=(const MclMatrix &Matrix)
236     {
237         return *this=*this+Matrix;
238     }
239     MclMatrix operator-=(const MclMatrix &Matrix)
240     {
241         return *this=*this-Matrix;
242     }
243 
244     MclMatrix operator*(double x) const
245     {
246         MclMatrix ans=*this;
247         for(int i=1;i<=row;i++)
248         {
249             for(int j=1;j<=col;j++)
250             {
251                 ans[i][j]*=x;
252             }
253         }
254         return ans;
255     }
256 
257 };
258 
259 MclMatrix vectorMulVector(const MclVector &A,const MclVector& B)
260 {
261     if(A.type==0)
262     {
263         MclMatrix ans=MclMatrix(A.n,B.n);
264         for(int i=1;i<=A.n;i++)
265         {
266             for(int j=1;j<=B.n;j++)
267             {
268                 ans[i][j]+=A[i]*B[j];
269             }
270         }
271         return ans;
272     }
273     else
274     {
275         assert(A.n==B.n);
276         MclMatrix ans=MclMatrix(1,1);
277         for(int i=1;i<=A.n;i++)
278         {
279             ans[1][1]+=A[i]*B[i];
280         }
281         return ans;
282     }
283 }
284 
285 int sgn(double x)
286 {
287     if(x<-0.000001) return -1;
288     if(x>0.000001) return 1;
289     return 0;
290 }
291 
292 /**
293   将矩阵A分解为一个正交矩阵Q和一个上三角矩阵R
294   A为任意实数矩阵
295 **/
296 std::pair<MclMatrix,MclMatrix> QRSplit(const MclMatrix &A)
297 {
298     assert(A.col==A.row);
299     int n=A.row;
300     MclMatrix Q=MclMatrix(n,n); Q.setIdentityMatrix();
301     MclMatrix R=A;
302 
303 
304     for(int i=1;i<n;i++)
305     {
306         MclVector s=R.GetCol(i);
307 
308         for(int j=1;j<i;j++) s[j]=0;
309 
310         if(sgn(s.length())==0) continue;
311 
312 
313         double c=s.length();
314         if(sgn(R[i][i])!=0) c*=-sgn(R[i][i]);
315         MclVector u=s; u[i]-=c;
316         MclVector uT=s.GetTranspose();
317 
318         MclMatrix H=MclMatrix(n,n);
319         H.setIdentityMatrix();
320 
321         H=H-vectorMulVector(u,uT)*(2.0/(u.length()*u.length()));
322         R=H*R;
323         Q=Q*H;
324     }
325     return std::make_pair(Q,R);
326 }