实对称矩阵不一定可逆。
如果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
Sw实对称,转为可逆:SW=SW+0.0001*trace(SW)*eye()
inv(Sw)不是对称矩阵,所以inv(Sw)*Sb不是实对称矩阵,所以opencv中不能直接用eigen()函数(只能处理对称矩阵)分解,而是用LDA.cpp中的一个特征分解类。
-----matlab可以直接用eigen()分解。
如果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