实对称矩阵的一些小知识——在LDA和LPP中遇到的

来源:互联网 发布:安藤政信 知乎 编辑:程序博客网 时间:2024/04/23 14:45
实对称矩阵不一定可逆。


Sw实对称,转为可逆:SW=SW+0.0001*trace(SW)*eye()

inv(Sw)不是对称矩阵,所以inv(Sw)*Sb不是实对称矩阵,所以opencv中不能直接用eigen()函数(只能处理对称矩阵)分解,而是用LDA.cpp中的一个特征分解类。

-----matlab可以直接用eigen()分解。


Sw和Sb都是实对称矩阵,是满秩的吗(特征值都不为0)?不一定!inv不一定可以用。
如果A是实矩阵,那么
1.AA^T是实对称半正定矩阵,AA^T特征值必定都大于等于0,但不一定满秩。

pinv(B)求的是矩阵B的Moore-Penrose逆,是B的一种广义逆,也就是你说的伪逆,该广义逆满足四个条件:
A*B*A = A
B*A*B = B
A*B 是海森矩阵
B*A是海森矩阵。
这个在矩阵论中有讲,你可以去看看



设A是对称矩阵
A^T = A
A^-1 = (A^T)^-1 = (A^-1)^T    (即A的逆也是对称矩阵)

可逆变换不改变秩

设n阶实对称矩阵A的秩为r(r<n),证明:存在可逆矩阵C,使得C^TAC=diag(d1,d2,...dr,0,...,0).其中di不等于0(

i=1,2,...)



 如果有n阶矩阵A,其各个元素都为实数,且aij=aji(转置为其本身),则称A为实对称矩阵。
  主要性质:
  1.实对称矩阵A的不同特征值所对应的特征向量是正交的。
  2.实对称矩阵A的特征值都是实数,特征向量都是实向量。
  3.n阶实对称矩阵A必可对角化。
  4.可用正交矩阵对角化。
  5.K重特征值必有K个线性无关的特征向量,或者说必有秩r(λE-A)=n-k。

对称:A的转置=A
正定:正惯性指数等于矩阵的阶数,所有特征值>0
而且正定矩阵肯定是对称的。B正定, 故inv(B)存在——B满秩



hermitian矩阵是实对称矩阵的推广,共轭转置等于本身的矩阵
A=A共轭转置
例如
 1  2i  3+i
-2i   5    6
 3-i   6    4



如果A是实矩阵,那么
1.AA^T是实对称半正定矩阵,AA^T特征值必定都大于0,但不一定满秩。
2.rank(A)=rank(AA^T)
3.||A||_2^2 = ||AA^T||_2
4.AA^T的特征值是A的奇异值的平方

只需要会证明4就行了,1,2,3都是推论。



里面有提到广义特征值的问题,但是书里面只介绍说(A+λM)a=0中的矩阵M是非奇异的,就是一般的矩阵。


使用eigs(A,B)或eig(A,B)时,B一定要正定(满秩是正定的条件之一)。且结果特征值正确,特征向量不一定对。
记住广义特征值一定可以用eig(inv(B)*A),且在LDA和LPP中更准确。

正定矩阵:
实对称、满秩、所有特征值大于0
(在大学线性代数教材范围内, 可认为正定矩阵都是对称矩阵

因为对正定矩阵的研究起源于对实二次型的研究, 矩阵是对应二次型的矩阵, 所以是对称的.
对复数域上的正定矩阵, 是共扼对称。之后又引入了广义正定矩阵, 且分有几类。)



*****************************
matlab的eig命令精度还是相当高的
在计算一般的矩阵,比如对称正定阵,其采用QR分解方法
在计算比较病态的或者比较恶心的矩阵比如非对称非正定其采用QZ分解方法,中间分解酉空间是采用迭代方法。
对于性质比较好的矩阵计算是很准的,但是对于比较病态的由于数值原因,无论什么算法,高阶特征值都是不很准的,

相对来说此时采用子空间迭代和laczos方法可能更好一些。
我个人编过各种算法计算特征值的程序(迭代和分解的都有),准确度上都没有matlab的eig高,所以总结下来matlab

的eig可以比较放心的使用

QR分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵,所以称为QR分解法,与此正规正交矩阵的通用符号Q有关。

MATLAB以qr函数来执行QR分解法, 其语法为[Q,R]=qr(A)。



广义特征值问题,即Ax=  Bx,

在Matlab中,使用eig()求解一般特征值问题和广义特征值。[V,D] = eig(A,B,flag), A和B时方阵,flag用来选择算

法,""qz""表示选择使用QZ算法。

也可以直接调用qz()来求解,[AA,BB,Q,Z,V] = qz(A,B,flag), flag 表示使用复数或实数计算,默认取值为复数。

在Lapack中,有四个函数都是用来求解广义特征值的,


?GEGS  Computes the generalized eigenvalues, Schur form, and left and/or right Schur vectors for a pair

of non-symmetric matrices.

?GGES  Computes the generalized eigenvalues, Schur form, and left and/or right Schur vectors for a pair

of non-symmetric matrices.

?GEGV  Computes the generalized eigenvalues, and left and/or right generalized eigenvectors for a pair

of non-symmetric matrices.

?GGEV  Computes the generalized eigenvalues, and left and/or right generalized eigenvectors for a pair

of non-symmetric matrices.



区别在于前两个分解之后会输出舒尔形式,后两个则输出广义特征向量。而且gegs和gegv都被gges和ggev代替。两个都

会用QZ分解求解广义特征值。

LAPACK也给出了QZ分解的函数dhgeqz,但要求输入H,T矩阵,对于一般的方阵,可以使用dgghrd将输入的方阵A,B变换

成H,T矩阵。
下面给出这四个函数的原型和测试程序。




matlab 矩阵函数
(2012-07-25 11:48:17)
转载▼
标签:
杂谈
    分类: MathWork
1.矩阵对角线元素的抽取
函数 diag
格式 X = diag(v,k)   %以向量v的元素作为矩阵X的第k条对角线元素,当k=0时,v为X的主对角线;当k>0时,v为上方

第k条对角线;当k<0时,v为下方第k条对角线。
X = diag(v)    %以v为主对角线元素,其余元素为0构成X。
v = diag(X,k)   %抽取X的第k条对角线元素构成向量v。k=0:抽取主对角线元素;k>0:抽取上方第k条对角线元素;

k<0抽取下方第k条对角线元素。
v = diag(X)    %抽取主对角线元素构成向量v。

2.上三角阵和下三角阵的抽取
函数 tril   %取下三角部分
格式 L = tril(X)     %抽取X的主对角线的下三角部分构成矩阵L
L = tril(X,k)    %抽取X的第k条对角线的下三角部分;k=0为主对角线;k>0为主对角线以上;k<0为主对角线以下。
函数 triu    %取上三角部分
格式 U = triu(X)    %抽取X的主对角线的上三角部分构成矩阵U
U = triu(X,k)   %抽取X的第k条对角线的上三角部分;k=0为主对角线;k>0为主对角线以上;k<0为主对角线以下。

3.矩阵的变维
矩阵的变维有两种方法,即用“:”和函数“reshape”,前者主要针对2个已知维数矩阵之间的变维操作;而后者是对

于一个矩阵的操作。
(1)“:”变维
(2)Reshape函数变维
格式 B = reshape(A,m,n)       %返回以矩阵A的元素构成的m×n矩阵B
B = reshape(A,m,n,p,…)   %将矩阵A变维为m×n×p×…
B = reshape(A,[m n p…])   %同上
B = reshape(A,siz)        %由siz决定变维的大小,元素个数与A中元素个数
相同。
(5)复制和平铺矩阵
函数 repmat
格式 B = repmat(A,m,n)       %将矩阵A复制m×n块,即B由m×n块A平铺而成。
B = repmat(A,[m n])      %与上面一致
B = repmat(A,[m n p…])   %B由m×n×p×…个A块平铺而成
repmat(A,m,n)           %当A是一个数a时,该命令产生一个全由a组成的m×n矩阵。

1.3 矩阵分解
1.3.1 Cholesky分解
函数 chol
格式 R = chol(X)      %如果X为n阶对称正定矩阵,则存在一个实的非奇异上三角阵R,满足R'*R = X;若X非正定,

则产生错误信息。
[R,p] = chol(X)   %不产生任何错误信息,若X为正定阵,则p=0,R与上相同;若X非正定,则p为正整数,R是有序的

上三角阵。
1.3.2 LU分解
矩阵的三角分解又称LU分解,它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积,即A=LU。
函数 lu
格式 [L,U] = lu(X)     %U为上三角阵,L为下三角阵或其变换形式,满足LU=X。
[L,U,P] = lu(X)   %U为上三角阵,L为下三角阵,P为单位矩阵的行变换矩阵,满足LU=PX。

1.3.3 QR分解
将矩阵A分解成一个正交矩阵与一个上三角矩阵的乘积。
函数 qr
格式 [Q,R] = qr(A)     %求得正交矩阵Q和上三角阵R,Q和R满足A=QR。
[Q,R,E] = qr(A)   %求得正交矩阵Q和上三角阵R,E为单位矩阵的变换形式,R的对角线元素按大小降序排列,满足

AE=QR。
[Q,R] = qr(A,0)    %产生矩阵A的“经济大小”分解
[Q,R,E] = qr(A,0) %E的作用是使得R的对角线元素降序,且Q*R=A(:, E)。
R = qr(A)         %稀疏矩阵A的分解,只产生一个上三角阵R,满足R'*R = A'*A,这种方法计算A'*A时减少了内在数

字信息的损耗。
[C,R] = qr(A,b)    %用于稀疏最小二乘问题:minimize||Ax-b||的两步解:[C,R] = qr(A,b),x = R\c。
R = qr(A,0)       %针对稀疏矩阵A的经济型分解
[C,R] = qr(A,b,0)   %针对稀疏最小二乘问题的经济型分解
函数 qrdelete
格式 [Q,R] = qrdelete(Q,R,j)   %返回将矩阵A的第j列移去后的新矩阵的qr分解
函数 qrinsert
格式 [Q,R] = qrinsert(Q,R,j,x)   %在矩阵A中第j列插入向量x后的新矩阵进行qr分解。若j大于A的列数,表示在A的

最后插入列x。
    
1.3.6 特征值分解
函数 eig
格式 d = eig(A)         %求矩阵A的特征值d,以向量形式存放d。
d = eig(A,B)       %A、B为方阵,求广义特征值d,以向量形式存放d。
[V,D] = eig(A)      %计算A的特征值对角阵D和特征向量V,使AV=VD成立。
[V,D] = eig(A,'nobalance')   %当矩阵A中有与截断误差数量级相差不远的值时,该指令可能更精确。'nobalance'起

误差调节作用。
[V,D] = eig(A,B)    %计算广义特征值向量阵V和广义特征值阵D,满足AV=BVD。
[V,D] = eig(A,B,flag)   % 由flag指定算法计算特征值D和特征向量V,flag的可能值为:'chol' 表示对B使用

Cholesky分解算法,这里A为对称Hermitian矩阵,B为正定阵。'qz' 表示使用QZ算法,这里A、B为非对称或非

Hermitian矩阵。
说明 一般特征值问题是求解方程: 解的问题。广义特征值问题是求方程: 解的问题。

*********************************

[a b]=eig(inv(B)*B)

a =

    0.8724   -0.0179   -0.8655
   -0.2186   -0.4471    0.4891
    0.4371    0.8943   -0.1084


b =

    1.0000         0         0
         0    1.0000         0
         0         0    1.0000


>> [a b]=eig(B*inv(B))

a =

    0.1294    0.6725   -0.3060
    0.9822    0.6101   -0.9485
    0.1364    0.4190    0.0824


b =

    1.0000         0         0
         0    1.0000         0
         0         0    1.0000


若B是满秩实对称矩阵,则BB-1和B-1B的特征值相同,特征向量不同。
B和3B的特征值倍增,特征向量相同。


若A是不满秩实对称矩阵,则,转为非奇异A=A+0.001*trace(A)*I或A=(A+A')/2,但结果有较大区别,不知道怎么用。

暂时A=A+0.001*trace(A)*I用于LDA和LPP中中Sw不满秩的问题,效果比A=(A+A')/2更佳。


d = eigs(A)     %求稀疏矩阵A的6个绝对值最大特征值d,d以向量形式存放。

d = eigs(A,B)       %求稀疏矩阵的广义特征值问题。满足AV=BVD,其中D为特征值对角阵,V为特征向量矩阵,B必须

是对称正定阵或Hermitian正定阵。

d = eigs(A,k)        %返回k个最大特征值

d = eigs(A,B,k)      %返回k个最大特征值





二 广义特征值
Ax = λBx 其中A,B均为n阶方阵。实际中往往 A 是Hermitian方阵(A=AH,不一定满秩),B是正定的Hermitian方阵

,A*V = B*V*D 是对应的广义特征值分解。

1) 一般的 V'AV = diag(不是由特征值组成的对角阵) , V'BV = diag(不是单位阵)。采用 inv(B)*Ax = λx

2)xi 和 xj 关于B矩阵带权正交时有如下性质: V'AV = D , V'BV = I(带权正交)。 采用第二种求法有此结论,优

先使用此方法

 MATLAB中是计算广义特征值[V,D] = eig(A,B) produces a diagonal matrix D of generalized eigenvalues and a

full matrix V whose columns are the corresponding eigenvectors so that A*V = B*V*D .

 A = [ 1 -1 1; -1 2 0; 1 0 3]; B = [5 2 -4;2 1 -2;-4 -2 5]; % B是 Hermitian 正定方阵

% 采用两种方法求广义特征值问题

% 1) B正定, 故inv(B)存在, 据Ax = λBx有 inv(B)*Ax = λx 将广义特征值问题转换为一般特征值问题求解λ和x

 [V, D] = eig(inv(B)*A)

V =

   -0.2280    0.8425    0.6300
    0.9525    0.4412    0.1543
    0.2019   -0.3090    0.7611
%即特征向量

D =

   23.0022         0         0
         0    0.0146         0
         0         0    2.9833

 [Q, R] = qr(V)

Q =

   -0.2280   -0.8927    0.3888
    0.9525   -0.2873   -0.1010
    0.2019    0.3473    0.9158


R =

    1.0000    0.1658    0.1570
         0   -0.9862   -0.3425
         0         0    0.9263
%%但是eig(A,B)计算结果却有问题
[V,D] = eig(A,B,'chol')

V =

    0.2918   -0.6339   -0.7162
    0.1528   -0.1553    2.9921
   -0.1070   -0.7658    0.6342


D =

    0.0146         0         0
         0    2.9833         0
         0         0   23.0022


% 2) B为Hermitian正定矩阵,故 B存在Cholesky分解 B = GG',其中G为下三角矩阵,将Ax = λBx变为两步求解λ和x

% step1:    inv(G)*A*(inv(G))' y = λ y 为一般特征值   
% step2: x = (inv(G))' y

 G = chol(B,'lower')

G =

    2.2361         0         0
    0.8944    0.4472         0
   -1.7889   -0.8944    1.0000

 S = inv(G)*A*inv(G)'

 [V, D] = eig(S)

V =

    0.0598   -0.9806    0.1865
   -0.7709   -0.1641   -0.6155
   -0.6342    0.1070    0.7658


D =

   23.0022         0         0
         0    0.0146         0
         0         0    2.9833

V_temp = inv(G)'*V   %V_temp中的各个列是 Ax = λBx 的特征向量

V_temp =

    0.7162   -0.2918    0.6339
   -2.9921   -0.1528    0.1553
   -0.6342    0.1070    0.7658

% 验证 V'AV = D , V'BV = I(带权正交) , 此时V中的各个广义特征向量有带权正交性

 V_temp'*B*V_temp

ans =

    1.0000   -0.0000   -0.0000
   -0.0000    1.0000   -0.0000
   -0.0000    0.0000    1.0000

V_temp'*A*V_temp  % 看看不就是广义特征值对角矩阵吗?

ans =

   23.0022   -0.0000    0.0000
   -0.0000    0.0146    0.0000
    0.0000    0.0000    2.9833

 

[Q, R] = qr(inv(G)'*V)

Q =

   -0.2280   -0.8927    0.3888
    0.9525   -0.2873   -0.1010
    0.2019    0.3473    0.9158


R =

   -3.1413   -0.0574    0.1580
         0    0.3416   -0.3446
         0         0    0.9320