特征值分解、奇異值分解、PCA概念整理

標簽： PCA特征值及向量奇異值QR算法 2014-01-18 14:21 16402人閱讀 評論(5) 收藏 舉報 本文章已收錄于： 分類： 機器學習（9） 作者同類文章X

版權聲明：本文為博主原創文章，未經博主允許不得轉載。

本文將分別介紹特征值分解、奇異值分解、及PCA的相關理論概念。

文章末尾將給出Householder矩陣變換、QR算法求解特征值、特征向量的代碼

其中，特征值分解、奇異值分解的相關內容，轉載自：

http://www.cnblogs.com/LeftNotEasy/archive/2011/01/19/svd-and-applications.html

考慮到本文50%以上的部分不是那個哥們的博客原文，所以我搞成原創標題了。。。。

一、特征值與特征向量的幾何意義

1.?????矩陣乘法

在介紹特征值與特征向量的幾何意義之前，先介紹矩陣乘法的幾何意義。

矩陣乘法對應了一個變換，是把任意一個向量變成另一個方向或長度的新向量。在這個變化過程中，原向量主要發生旋轉、伸縮的變化。如果矩陣對某些向量只發生伸縮變換，不產生旋轉效果，那么這些向量就稱為這個矩陣的特征向量，伸縮的比例就是特征值。

比如：，它對應的線性變換是下面的形式形式：

因為，這個矩陣乘以一個向量(x,y)的結果是：。由于矩陣M是對稱的，所以這個變換是一個對 x , y 軸的一個拉伸變換。【當M中元素值大于1時，是拉伸；當值小于1時，是縮短】

那么如果矩陣M不是對稱的，比如：，它所描述的變換如下圖所示：

這其實是在平面上對一個軸進行的拉伸變換【如藍色箭頭所示】，在圖中藍色箭頭是一個最主要的變化方向。變化方向可能有不止一個，但如果我們想要描述好一個變換，那我們就描述好這個變換主要的變化方向就好了。

2.?????特征值分解與特征向量

如果說一個向量v是方陣A的特征向量，將一定可以表示成下面的形式：

λ為特征向量 v 對應的特征值。特征值分解是將一個矩陣分解為如下形式：

其中，Q是這個矩陣A的特征向量組成的矩陣，Σ是一個對角矩陣，每一個對角線元素就是一個特征值，里面的特征值是由大到小排列的，這些特征值所對應的特征向量就是描述這個矩陣變化方向（從主要的變化到次要的變化排列）。也就是說矩陣A的信息可以由其特征值和特征向量表示。

對于矩陣為高維的情況下，那么這個矩陣就是高維空間下的一個線性變換。可以想象，這個變換也同樣有很多的變換方向，我們通過特征值分解得到的前N個特征向量，那么就對應了這個矩陣最主要的N個變化方向。我們利用這前N個變化方向，就可以近似這個矩陣（變換）。

總結一下，特征值分解可以得到特征值與特征向量，特征值表示的是這個特征到底有多重要，而特征向量表示這個特征是什么。不過，特征值分解也有很多的局限，比如說變換的矩陣必須是方陣。

二、奇異值分解

1.?????奇異值

特征值分解是一個提取矩陣特征很不錯的方法，但是它只是對方陣而言的，在現實的世界中，我們看到的大部分矩陣都不是方陣，比如說有N個學生，每個學生有M科成績，這樣形成的一個N * M的矩陣就不可能是方陣，我們怎樣才能描述這樣普通的矩陣呢的重要特征呢？奇異值分解可以用來干這個事情，奇異值分解是一個能適用于任意的矩陣的一種分解的方法：

分解形式：

假設A是一個M * N的矩陣，那么得到的U是一個M * M的方陣（稱為左奇異向量），Σ是一個M * N的矩陣（除了對角線的元素都是0，對角線上的元素稱為奇異值），V^T(V的轉置)是一個N * N的矩陣（稱為右奇異向量）。

2.?????奇異值與特征值

那么奇異值和特征值是怎么對應起來的呢？我們將一個矩陣A的轉置乘以 A，并對(AA^T)求特征值，則有下面的形式：

這里V就是上面的右奇異向量，另外還有：

這里的σ就是奇異值，u就是上面說的左奇異向量。【證明那個哥們也沒給】

奇異值σ跟特征值類似，在矩陣Σ中也是從大到小排列，而且σ的減少特別的快，在很多情況下，前10%甚至1%的奇異值的和就占了全部的奇異值之和的99%以上了。也就是說，我們也可以用前r（ r遠小于m、n ）個的奇異值來近似描述矩陣，即部分奇異值分解：

右邊的三個矩陣相乘的結果將會是一個接近于A的矩陣，在這兒，r越接近于n，則相乘的結果越接近于A。

三、PCA主成份分析

主成分分析(PrincipalComponents Analysis。即PCA，也稱為K-L變換)，是圖像壓縮中的一種最優正交變換。PCA用于統計特征提取構成了子空間法模式識別的基礎。它從圖像整體代數特征出發，基于圖像的總體信息進行分類識別。PCA的核心思想是利用較少數量的特征對樣本進行描述以達到降低特征空間維數的目的。

1.? PCA理論

給定一副N*N大小圖像，將它表示成一個N²*1維向量，向量中元素為像素點灰度，按行存儲，如下列公式分別表示第i張圖片和n張圖片平均值：

令N²*n矩陣X為：

注意，矩陣減去平均值相當于將坐標系原點移動到平均值位置。

設Q=XX^T，則Q是一個N²* N²矩陣：

,Q是方陣

,Q是對稱矩陣。

,Q被成為協方差矩陣，

,Q的數據量非常龐大

??? 那么，X中的每個元素x_j可以被如下表達：

其中,e_i是Q中非零特征值對應的特征向量。由特征向量e₁,e₂,…,e_n組成的空間叫做張成特征空間。對于N*N圖像，e₁,e₂,…,e_n是N²*1維相互正交的向量。尺度g_ji是x_j在空間中的坐標。

2.? 實現PCA

為了降維，可以對特征值設定閾值或按照其他準則，尋找協方差矩陣Q中前k個特征向量。這里Q十分龐大，對于一副256*256的圖像，Q的大小為65536*65536！替代方案是，考慮矩陣

.P和Q都是對稱矩陣

.P≠Q^T

.Q的大小是N²*N²，而P大小為n*n

.n為訓練樣本圖像數量，通常n<<N

設e是矩陣P的特征值λ對應的特征向量，則有：

這里，X*e也是矩陣Q的特征值λ對應的特征向量【這是用求特征值分解方法，下面介紹用SVD方法】

3.? PCA與奇異值分解SVD

任何一個m*n矩陣都能進行奇異值分解，拆分為3個矩陣相乘的形式。由于SVD得出的奇異向量也是從奇異值由大到小排列的，按PCA的觀點來看，就是方差最大的坐標軸就是第一個奇異向量，方差次大的坐標軸就是第二個奇異向量…。我們可以對Q進行奇異值分解。

.U就是QQ^T的特征向量

.V就是Q^TQ的特征向量

.D中奇異值的平方就是QQ^T和Q^TQ的特征值

=======================================================================================================================

上面講了一大堆，就是為了下一篇PCA人臉識別做鋪墊的，給你一副圖像，要從圖像庫中得到匹配的圖像，怎么弄？如果是兩兩做像素點比較是不可能完成的任務，時間上廢掉了。如果用其他特征點代替也許可以，但容易漏檢吧，這邊不擴展。我們必須對圖像數據的協方差矩陣進行降維，所以用到了PCA。

而具體如何實現PCA呢？關鍵是特征值及相應特征向量的求取。matlab有個eig函數，OpenCV也有相應的函數。由于不想被別人牽制，我自己查了資料，發現QR算法可以用來求實對稱矩陣的全部特征值和特征向量。【雅可比算法也可以，就是速度太慢了；而上面介紹的SVD實現PCA還沒見過，文獻上說SVD和PCA是等價的】

=======================================================================================================================

以下內容，來自《C常用算法程序集第二版》，這絕對是搞科研的好書！

在用QR算法求解特征值和向量之前，必須將實對稱矩陣轉化為三對角矩陣。【由于我們的協方差矩陣是實對稱矩陣，因此不用轉化為Hessen berg矩陣，QR算法是一個迭代的過程，具體算法太長了，我不貼出來了，有需要的，自己去下載這本書的PDF文檔或其他資料】

1.約化對稱矩陣為三對角矩陣的Householder變換法：

例：

【其他高維矩陣也行，大家可以把數據存在txt文本中，然后讀取進來】

代碼：

[cpp] view plaincopy print?

//?HouseHolder_Transform.cpp?:?定義控制臺應用程序的入口點。??

//??

??

#include?"stdafx.h"??

#include?"math.h"??

??

void?cstrq(double?a[],int?n,double?q[],double?b[],double?c[]);??

??

int?_tmain(int?argc,?_TCHAR*?argv[])??

{??

????int?i,j;??

????static?double?b[5],c[5],q[25];??

????static?double?a[25]?=?{10.0,1.0,2.0,3.0,4.0,1.0,9.0,-1.0,2.0,-3.0,2.0,-1.0,7.0,3.0,-5.0,3.0,2.0,3.0,12.0,-1.0,4.0,-3.0,-5.0,-1.0,15.0};??

????cstrq(a,5,q,b,c);??

??

????printf("MAT?A?is:\n");??

????for?(i=0;i<5;i++)??

????{??

????????for?(j=0;j<5;j++)??

????????{??

????????????printf("%13.7e?",a[i*5+j]);??

????????}??

????????printf("\n");??

????}??

????printf("\n");??

????printf("MAT?Q?is:\n");??

????for?(i=0;i<5;i++)??

????{??

????????for?(j=0;j<5;j++)??

????????{??

????????????printf("%13.7e?",q[i*5+j]);??

????????}??

????????printf("\n");??

????}??

????printf("\n");??

????printf("MAT?B?is:\n");??

????for?(i=0;i<5;i++)??

????{??

????????printf("%13.7e?",b[i]);??

????}??

????printf("\n\n");??

????printf("MAT?C?is:\n");??

????for?(i=0;i<5;i++)??

????{??

????????printf("%13.7e?",c[i]);??

????}??

????printf("\n\n");??

????return?0;??

}??

??

??

??

void?cstrq(double?a[],int?n,double?q[],double?b[],double?c[])??

{??

????int?i,j,k,u,v;??

????double?h,f,g,h2;??

????for?(i=0;?i<=n-1;?i++)??

????????for?(j=0;?j<=n-1;?j++)??

????????{?u=i*n+j;?q[u]=a[u];}??

????????for?(i=n-1;?i>=1;?i--)??

????????{?h=0.0;??

????????if?(i>1)??

????????????for?(k=0;?k<=i-1;?k++)??

????????????{?u=i*n+k;?h=h+q[u]*q[u];}??

????????????if?(h+1.0==1.0)??

????????????{?c[i]=0.0;??

????????????if?(i==1)?c[i]=q[i*n+i-1];??

????????????b[i]=0.0;??

????????????}??

????????????else??

????????????{?c[i]=sqrt(h);??

????????????u=i*n+i-1;??

????????????if?(q[u]>0.0)?c[i]=-c[i];??

????????????h=h-q[u]*c[i];??

????????????q[u]=q[u]-c[i];??

????????????f=0.0;??

????????????for?(j=0;?j<=i-1;?j++)??

????????????{?q[j*n+i]=q[i*n+j]/h;??

????????????g=0.0;??

????????????for?(k=0;?k<=j;?k++)??

????????????????g=g+q[j*n+k]*q[i*n+k];??

????????????if?(j+1<=i-1)??

????????????????for?(k=j+1;?k<=i-1;?k++)??

????????????????????g=g+q[k*n+j]*q[i*n+k];??

????????????c[j]=g/h;??

????????????f=f+g*q[j*n+i];??

????????????}??

????????????h2=f/(h+h);??

????????????for?(j=0;?j<=i-1;?j++)??

????????????{?f=q[i*n+j];??

????????????g=c[j]-h2*f;??

????????????c[j]=g;??

????????????for?(k=0;?k<=j;?k++)??

????????????{?u=j*n+k;??

????????????q[u]=q[u]-f*c[k]-g*q[i*n+k];??

????????????}??

????????????}??

????????????b[i]=h;??

????????????}??

????????}??

????????for?(i=0;?i<=n-2;?i++)?c[i]=c[i+1];??

????????c[n-1]=0.0;??

????????b[0]=0.0;??

????????for?(i=0;?i<=n-1;?i++)??

????????{?if?((b[i]!=0.0)&&(i-1>=0))??

????????for?(j=0;?j<=i-1;?j++)??

????????{?g=0.0;??

????????for?(k=0;?k<=i-1;?k++)??

????????????g=g+q[i*n+k]*q[k*n+j];??

????????for?(k=0;?k<=i-1;?k++)??

????????{?u=k*n+j;??

????????q[u]=q[u]-g*q[k*n+i];??

????????}??

????????}??

????????u=i*n+i;??

????????b[i]=q[u];?q[u]=1.0;??

????????if?(i-1>=0)??

????????????for?(j=0;?j<=i-1;?j++)??

????????????{?q[i*n+j]=0.0;?q[j*n+i]=0.0;}??

????????}??

????????return;??

}??

// HouseHolder_Transform.cpp : 定義控制臺應用程序的入口點。 //#include "stdafx.h" #include "math.h"void cstrq(double a[],int n,double q[],double b[],double c[]);int _tmain(int argc, _TCHAR* argv[]) {int i,j;static double b[5],c[5],q[25];static double a[25] = {10.0,1.0,2.0,3.0,4.0,1.0,9.0,-1.0,2.0,-3.0,2.0,-1.0,7.0,3.0,-5.0,3.0,2.0,3.0,12.0,-1.0,4.0,-3.0,-5.0,-1.0,15.0};cstrq(a,5,q,b,c);printf("MAT A is:\n");for (i=0;i<5;i++){for (j=0;j<5;j++){printf("%13.7e ",a[i*5+j]);}printf("\n");}printf("\n");printf("MAT Q is:\n");for (i=0;i<5;i++){for (j=0;j<5;j++){printf("%13.7e ",q[i*5+j]);}printf("\n");}printf("\n");printf("MAT B is:\n");for (i=0;i<5;i++){printf("%13.7e ",b[i]);}printf("\n\n");printf("MAT C is:\n");for (i=0;i<5;i++){printf("%13.7e ",c[i]);}printf("\n\n");return 0; }void cstrq(double a[],int n,double q[],double b[],double c[]) {int i,j,k,u,v;double h,f,g,h2;for (i=0; i<=n-1; i++)for (j=0; j<=n-1; j++){ u=i*n+j; q[u]=a[u];}for (i=n-1; i>=1; i--){ h=0.0;if (i>1)for (k=0; k<=i-1; k++){ u=i*n+k; h=h+q[u]*q[u];}if (h+1.0==1.0){ c[i]=0.0;if (i==1) c[i]=q[i*n+i-1];b[i]=0.0;}else{ c[i]=sqrt(h);u=i*n+i-1;if (q[u]>0.0) c[i]=-c[i];h=h-q[u]*c[i];q[u]=q[u]-c[i];f=0.0;for (j=0; j<=i-1; j++){ q[j*n+i]=q[i*n+j]/h;g=0.0;for (k=0; k<=j; k++)g=g+q[j*n+k]*q[i*n+k];if (j+1<=i-1)for (k=j+1; k<=i-1; k++)g=g+q[k*n+j]*q[i*n+k];c[j]=g/h;f=f+g*q[j*n+i];}h2=f/(h+h);for (j=0; j<=i-1; j++){ f=q[i*n+j];g=c[j]-h2*f;c[j]=g;for (k=0; k<=j; k++){ u=j*n+k;q[u]=q[u]-f*c[k]-g*q[i*n+k];}}b[i]=h;}}for (i=0; i<=n-2; i++) c[i]=c[i+1];c[n-1]=0.0;b[0]=0.0;for (i=0; i<=n-1; i++){ if ((b[i]!=0.0)&&(i-1>=0))for (j=0; j<=i-1; j++){ g=0.0;for (k=0; k<=i-1; k++)g=g+q[i*n+k]*q[k*n+j];for (k=0; k<=i-1; k++){ u=k*n+j;q[u]=q[u]-g*q[k*n+i];}}u=i*n+i;b[i]=q[u]; q[u]=1.0;if (i-1>=0)for (j=0; j<=i-1; j++){ q[i*n+j]=0.0; q[j*n+i]=0.0;}}return; }
計算結果：

即上述計算結果返回的三對角陣T為：

2.下面，我們將在三對角矩陣的基礎上使用QR算法計算全部特征值和特征向量

例，同樣對上面那個5階矩陣，先求三對角矩陣，再求其全部特征值和特征向量

最大迭代次數為60，誤差為0.000001

代碼：

[cpp] view plaincopy print?

#include?"stdafx.h"??

#include?"math.h"??

??

void?cstrq(double?a[],int?n,double?q[],double?b[],double?c[]);??

int?csstq(int?n,double?b[],double?c[],double?q[],double?eps,int?l);??

??

int?_tmain(int?argc,?_TCHAR*?argv[])??

{??

????int?i,j,k,l=60;??

????static?double?b[5],c[5],q[25];??

????static?double?a[25]?=?{10.0,1.0,2.0,3.0,4.0,1.0,9.0,-1.0,2.0,-3.0,2.0,-1.0,7.0,3.0,-5.0,3.0,2.0,3.0,12.0,-1.0,4.0,-3.0,-5.0,-1.0,15.0};??

????double?eps?=?0.000001;??

??

????cstrq(a,5,q,b,c);??

????k?=?csstq(5,b,c,q,eps,l);??

??

????printf("MAT?A?is:\n");??

????for?(i=0;i<5;i++)??

????{??

????????for?(j=0;j<5;j++)??

????????{??

????????????printf("%13.7e?",a[i*5+j]);??

????????}??

????????printf("\n");??

????}??

????printf("\n");??

????printf("MAT?B?is:\n");??

????for?(i=0;i<5;i++)??

????{??

????????printf("%13.7e?",b[i]);??

????}??

????printf("\n\n");??

????printf("MAT?Q?is:\n");??

????for?(i=0;i<5;i++)??

????{??

????????for?(j=0;j<5;j++)??

????????{??

????????????printf("%13.7e?",q[i*5+j]);??

????????}??

????????printf("\n");??

????}??

????printf("\n");??

????return?0;??

}??

??

??

??

void?cstrq(double?a[],int?n,double?q[],double?b[],double?c[])??

{??

????int?i,j,k,u,v;??

????double?h,f,g,h2;??

????for?(i=0;?i<=n-1;?i++)??

????????for?(j=0;?j<=n-1;?j++)??

????????{?u=i*n+j;?q[u]=a[u];}??

????????for?(i=n-1;?i>=1;?i--)??

????????{?h=0.0;??

????????if?(i>1)??

????????????for?(k=0;?k<=i-1;?k++)??

????????????{?u=i*n+k;?h=h+q[u]*q[u];}??

????????????if?(h+1.0==1.0)??

????????????{?c[i]=0.0;??

????????????if?(i==1)?c[i]=q[i*n+i-1];??

????????????b[i]=0.0;??

????????????}??

????????????else??

????????????{?c[i]=sqrt(h);??

????????????u=i*n+i-1;??

????????????if?(q[u]>0.0)?c[i]=-c[i];??

????????????h=h-q[u]*c[i];??

????????????q[u]=q[u]-c[i];??

????????????f=0.0;??

????????????for?(j=0;?j<=i-1;?j++)??

????????????{?q[j*n+i]=q[i*n+j]/h;??

????????????g=0.0;??

????????????for?(k=0;?k<=j;?k++)??

????????????????g=g+q[j*n+k]*q[i*n+k];??

????????????if?(j+1<=i-1)??

????????????????for?(k=j+1;?k<=i-1;?k++)??

????????????????????g=g+q[k*n+j]*q[i*n+k];??

????????????c[j]=g/h;??

????????????f=f+g*q[j*n+i];??

????????????}??

????????????h2=f/(h+h);??

????????????for?(j=0;?j<=i-1;?j++)??

????????????{?f=q[i*n+j];??

????????????g=c[j]-h2*f;??

????????????c[j]=g;??

????????????for?(k=0;?k<=j;?k++)??

????????????{?u=j*n+k;??

????????????q[u]=q[u]-f*c[k]-g*q[i*n+k];??

????????????}??

????????????}??

????????????b[i]=h;??

????????????}??

????????}??

????????for?(i=0;?i<=n-2;?i++)?c[i]=c[i+1];??

????????c[n-1]=0.0;??

????????b[0]=0.0;??

????????for?(i=0;?i<=n-1;?i++)??

????????{?if?((b[i]!=0.0)&&(i-1>=0))??

????????for?(j=0;?j<=i-1;?j++)??

????????{?g=0.0;??

????????for?(k=0;?k<=i-1;?k++)??

????????????g=g+q[i*n+k]*q[k*n+j];??

????????for?(k=0;?k<=i-1;?k++)??

????????{?u=k*n+j;??

????????q[u]=q[u]-g*q[k*n+i];??

????????}??

????????}??

????????u=i*n+i;??

????????b[i]=q[u];?q[u]=1.0;??

????????if?(i-1>=0)??

????????????for?(j=0;?j<=i-1;?j++)??

????????????{?q[i*n+j]=0.0;?q[j*n+i]=0.0;}??

????????}??

????????return;??

}??

??

int?csstq(int?n,double?b[],double?c[],double?q[],double?eps,int?l)??

{??

????int?i,j,k,m,it,u,v;??

????double?d,f,h,g,p,r,e,s;??

????c[n-1]=0.0;?d=0.0;?f=0.0;??

????for?(j=0;?j<=n-1;?j++)??

????{?it=0;??

????h=eps*(fabs(b[j])+fabs(c[j]));??

????if?(h>d)?d=h;??

????m=j;??

????while?((m<=n-1)&&(fabs(c[m])>d))?m=m+1;??

????if?(m!=j)??

????{?do??

????{?if?(it==l)??

????{?printf("fail\n");??

????return(-1);??

????}??

????it=it+1;??

????g=b[j];??

????p=(b[j+1]-g)/(2.0*c[j]);??

????r=sqrt(p*p+1.0);??

????if?(p>=0.0)?b[j]=c[j]/(p+r);??

????else?b[j]=c[j]/(p-r);??

????h=g-b[j];??

????for?(i=j+1;?i<=n-1;?i++)??

????????b[i]=b[i]-h;??

????f=f+h;?p=b[m];?e=1.0;?s=0.0;??

????for?(i=m-1;?i>=j;?i--)??

????{?g=e*c[i];?h=e*p;??

????if?(fabs(p)>=fabs(c[i]))??

????{?e=c[i]/p;?r=sqrt(e*e+1.0);??

????c[i+1]=s*p*r;?s=e/r;?e=1.0/r;??

????}??

????else??

????{?e=p/c[i];?r=sqrt(e*e+1.0);??

????c[i+1]=s*c[i]*r;??

????s=1.0/r;?e=e/r;??

????}??

????p=e*b[i]-s*g;??

????b[i+1]=h+s*(e*g+s*b[i]);??

????for?(k=0;?k<=n-1;?k++)??

????{?u=k*n+i+1;?v=u-1;??

????h=q[u];?q[u]=s*q[v]+e*h;??

????q[v]=e*q[v]-s*h;??

????}??

????}??

????c[j]=s*p;?b[j]=e*p;??

????}??

????while?(fabs(c[j])>d);??

????}??

????b[j]=b[j]+f;??

????}??

????for?(i=0;?i<=n-1;?i++)??

????{?k=i;?p=b[i];??

????if?(i+1<=n-1)??

????{?j=i+1;??

????while?((j<=n-1)&&(b[j]<=p))??

????{?k=j;?p=b[j];?j=j+1;}??

????}??

????if?(k!=i)??

????{?b[k]=b[i];?b[i]=p;??

????for?(j=0;?j<=n-1;?j++)??

????{?u=j*n+i;?v=j*n+k;??

????p=q[u];?q[u]=q[v];?q[v]=p;??

????}??

????}??

????}??

????return(1);??

}??

#include "stdafx.h" #include "math.h"void cstrq(double a[],int n,double q[],double b[],double c[]); int csstq(int n,double b[],double c[],double q[],double eps,int l);int _tmain(int argc, _TCHAR* argv[]) {int i,j,k,l=60;static double b[5],c[5],q[25];static double a[25] = {10.0,1.0,2.0,3.0,4.0,1.0,9.0,-1.0,2.0,-3.0,2.0,-1.0,7.0,3.0,-5.0,3.0,2.0,3.0,12.0,-1.0,4.0,-3.0,-5.0,-1.0,15.0};double eps = 0.000001;cstrq(a,5,q,b,c);k = csstq(5,b,c,q,eps,l);printf("MAT A is:\n");for (i=0;i<5;i++){for (j=0;j<5;j++){printf("%13.7e ",a[i*5+j]);}printf("\n");}printf("\n");printf("MAT B is:\n");for (i=0;i<5;i++){printf("%13.7e ",b[i]);}printf("\n\n");printf("MAT Q is:\n");for (i=0;i<5;i++){for (j=0;j<5;j++){printf("%13.7e ",q[i*5+j]);}printf("\n");}printf("\n");return 0; }void cstrq(double a[],int n,double q[],double b[],double c[]) {int i,j,k,u,v;double h,f,g,h2;for (i=0; i<=n-1; i++)for (j=0; j<=n-1; j++){ u=i*n+j; q[u]=a[u];}for (i=n-1; i>=1; i--){ h=0.0;if (i>1)for (k=0; k<=i-1; k++){ u=i*n+k; h=h+q[u]*q[u];}if (h+1.0==1.0){ c[i]=0.0;if (i==1) c[i]=q[i*n+i-1];b[i]=0.0;}else{ c[i]=sqrt(h);u=i*n+i-1;if (q[u]>0.0) c[i]=-c[i];h=h-q[u]*c[i];q[u]=q[u]-c[i];f=0.0;for (j=0; j<=i-1; j++){ q[j*n+i]=q[i*n+j]/h;g=0.0;for (k=0; k<=j; k++)g=g+q[j*n+k]*q[i*n+k];if (j+1<=i-1)for (k=j+1; k<=i-1; k++)g=g+q[k*n+j]*q[i*n+k];c[j]=g/h;f=f+g*q[j*n+i];}h2=f/(h+h);for (j=0; j<=i-1; j++){ f=q[i*n+j];g=c[j]-h2*f;c[j]=g;for (k=0; k<=j; k++){ u=j*n+k;q[u]=q[u]-f*c[k]-g*q[i*n+k];}}b[i]=h;}}for (i=0; i<=n-2; i++) c[i]=c[i+1];c[n-1]=0.0;b[0]=0.0;for (i=0; i<=n-1; i++){ if ((b[i]!=0.0)&&(i-1>=0))for (j=0; j<=i-1; j++){ g=0.0;for (k=0; k<=i-1; k++)g=g+q[i*n+k]*q[k*n+j];for (k=0; k<=i-1; k++){ u=k*n+j;q[u]=q[u]-g*q[k*n+i];}}u=i*n+i;b[i]=q[u]; q[u]=1.0;if (i-1>=0)for (j=0; j<=i-1; j++){ q[i*n+j]=0.0; q[j*n+i]=0.0;}}return; }int csstq(int n,double b[],double c[],double q[],double eps,int l) {int i,j,k,m,it,u,v;double d,f,h,g,p,r,e,s;c[n-1]=0.0; d=0.0; f=0.0;for (j=0; j<=n-1; j++){ it=0;h=eps*(fabs(b[j])+fabs(c[j]));if (h>d) d=h;m=j;while ((m<=n-1)&&(fabs(c[m])>d)) m=m+1;if (m!=j){ do{ if (it==l){ printf("fail\n");return(-1);}it=it+1;g=b[j];p=(b[j+1]-g)/(2.0*c[j]);r=sqrt(p*p+1.0);if (p>=0.0) b[j]=c[j]/(p+r);else b[j]=c[j]/(p-r);h=g-b[j];for (i=j+1; i<=n-1; i++)b[i]=b[i]-h;f=f+h; p=b[m]; e=1.0; s=0.0;for (i=m-1; i>=j; i--){ g=e*c[i]; h=e*p;if (fabs(p)>=fabs(c[i])){ e=c[i]/p; r=sqrt(e*e+1.0);c[i+1]=s*p*r; s=e/r; e=1.0/r;}else{ e=p/c[i]; r=sqrt(e*e+1.0);c[i+1]=s*c[i]*r;s=1.0/r; e=e/r;}p=e*b[i]-s*g;b[i+1]=h+s*(e*g+s*b[i]);for (k=0; k<=n-1; k++){ u=k*n+i+1; v=u-1;h=q[u]; q[u]=s*q[v]+e*h;q[v]=e*q[v]-s*h;}}c[j]=s*p; b[j]=e*p;}while (fabs(c[j])>d);}b[j]=b[j]+f;}for (i=0; i<=n-1; i++){ k=i; p=b[i];if (i+1<=n-1){ j=i+1;while ((j<=n-1)&&(b[j]<=p)){ k=j; p=b[j]; j=j+1;}}if (k!=i){ b[k]=b[i]; b[i]=p;for (j=0; j<=n-1; j++){ u=j*n+i; v=j*n+k;p=q[u]; q[u]=q[v]; q[v]=p;}}}return(1); }
這里，又把householder貼了一遍。。。

計算結果：

這里，我們要注意：

數組q中第j列為數組b中第j個特征值對應的特征向量

頂

5

踩

0

• 上一篇Delaunay 三角網格學習
• 下一篇PCA人臉識別學習及C語言實現

我的同類文章

機器學習（9） http://blog.csdn.net

• ?Gauss-Newton算法學習2016-06-08閱讀2927
• ?機器學習之旅---奇異值分解2014-11-22閱讀2946
• ?機器學習之旅---logistic回歸2014-10-12閱讀4802
• ?離散隨機線性系統的卡爾曼濾波器基本原理及實現2014-06-13閱讀5079
• ?SVM理論部分介紹2013-12-28閱讀2716

• ?《Master Opencv...讀書筆記》非剛性人臉跟蹤 III2015-02-28閱讀3040
• ?機器學習之旅---SVM分類器2014-11-07閱讀27607
• ?機器學習之旅---樸素貝葉斯分類器2014-09-24閱讀2372
• ?matlab體驗svm算法【非實現】2013-12-30閱讀4855

總結

以上是生活随笔為你收集整理的特征值 奇异值分解 概念整理的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。