矩阵运算

1. 矩阵运算

矩阵运算包括矩阵分析、特征值求解、奇异值等

1.1 矩阵分析

1.1.0 向量和矩阵的范数运算

求向量范数的函数

N=norm(x,p)：对任意大于1的p值，返回向量x的p阶范数
N=norm(x)：返回向量的2阶范数，相当于N=norm(x,2)
N=norm(x,inf)：返回向量的$\infty$阶范数，相当于N=max(abs(x))
N=norm(x,-inf)：返回向量的–$\infty$阶范数，相当于N=min(abs(x))

求矩阵范数的函数

N=norm(A)：计算矩阵的2阶范数，也就是最大奇异值
N=norm(A,p)：根据参数p的值不同，求不同阶的范数值
当p=1时，计算矩阵A的1阶范数，相当于max(sum(abs(A)))
当p=2时，计算矩阵A的2阶范数，相当于norm(A)
当p=inf时，计算矩阵A的$\infty$阶范数，相当于max(sum(abs(A’)))
当p=pro时，计算矩阵A的F范数(Frobenius范数)，相当于sqrt(sum(diag(A’*A)))

求向量x的2阶范数范例：

当矩阵维数比较大时，会导致计算矩阵范数的时间比较长，并且当一个近似的范数值满足要求时，可以考虑用函数 normest( )来估计2阶范数值。
函数 normest( ) 最初是为了提供给稀疏矩阵(sparse matrix)使用的，同时它也能接受满矩阵的输入，一般在满矩阵维数比较大时使用

normest(S)：估计矩阵S的2阶范数值
normest(S,tol)：使用 tol 作为允许的相对误差

求矩阵的范式示例：

1.1.1 矩阵的秩

矩阵A中线性无关的列向量个数成为列秩，线性无关的行向量个数称为行秩

rank(A)：默认允许误差计算矩阵的秩
rank(A,tol)：给定允许误差计算矩阵的秩，tol=max(size(A))*eps(norm(A))

求矩阵的秩示例:

1.1.2 矩阵的行列式

det(A)：计算矩阵A的行列式

1.1.3 矩阵的迹

矩阵的迹定义为矩阵对角元素之和

trace(A)：矩阵Ａ的迹

1.1.4 矩阵的化零矩阵

Z=null(A)：返回矩阵A的一个化零矩阵，如果化零矩阵不存在则返回空矩阵
z=null(A,‘r’)：返回有理数形式的化零矩阵

求矩阵的化零矩阵范例：

1.1.5 矩阵的正交空间

矩阵A的正交空间Q具有Q’$\cdot$Q=I的性质，并且Q的列矢量构成的线性空间与矩阵A的列矢量构成的线性空间相同，且正交空间Q与矩阵A具有相同的秩

orth(A)：矩阵A的正交空间

1.1.6 矩阵的约化行阶梯形式

矩阵的约化行阶梯形式是高斯－约旦消去法解线性方程组的结果

R=rref(A)：返回矩阵A的约化行阶梯形式R
[R,jb]=rref(A)：返回矩阵A的约化行阶梯形式R，并返回1×r的向量jb
r为矩阵A的秩
A(:,jb)是矩阵A的列矢量构成的线性空间
R(1:r,jb)是r×r的单位矩阵
[R,jb]=rref(A,tol)：以tol作为允许的相对误差计算矩阵A的秩

求矩阵A的约化行阶梯形式示例：

1.1.7 矩阵空间之间的夹角

矩阵空间之间的夹角代表两个矩阵线性相关的程度
如果夹角很小，它们之间的线性相关度就很高
反之，它们之间的线性相关度不大

theta=subspace(A,B)：返回矩阵A和矩阵B之间的夹角

求矩阵A和B之间的夹角

1.2 矩阵分解

函数	功能描述
chol	Cholesky分解
cholinc	稀疏矩阵的不完全Cholesky分解
lu	矩阵LU分解｜
luinc	稀疏矩阵的不完全LU分解
qr	正交三角分解
svd	奇异值分解
gsvd	一般奇异值分解
schur	舒尔分解

线性方程组的求解主要基于4种基本的矩阵分解

1. Cholesky分解
2. 一般方阵的高斯消去法分解
3. 舒尔分解
4. 矩形矩阵的正交分解

1.2.0 对称正定矩阵的Cholesky分解

R=chol(X)：其中X为对称正定矩阵，R是上三角矩阵，使得X=R’$\cdot$R
[R,p]=chol(X)：
当X是正定矩阵时，返回的上三角矩阵R满足X=R’$\cdot$R，且p=0；
当X是非正定矩阵时，返回值p是正整数，R是上三角矩阵，其阶数p-1，且满足X(1:p-1,1:p-1)=R’$\cdot$R

考虑线性方程组$Ax=b$，其中A可以做Cholesky分解，使得$A=R^{\prime }\cdot R$，这样线性方程组就可以改写成$R^{\prime }\cdot R\cdot x=b$，由于左除算符( \ )可以快速处理三角矩阵，因此得出：$$x=R\backslash (R^{\prime }\backslash b)$$如果A是方阵$n\times n$的方阵，则chol(Ａ)的计算复杂度 $O(n^3)$，而左除算符的计算复杂度只有 $O(n^2)$

利用chol进行矩阵分解范例：

在矩阵中，若数值为0的元素数目远远多于非0元素的数目，并且非0元素分布没有规律时，则称该矩阵为稀疏矩阵；
与之相反，若非0元素数目占大多数时，则称该矩阵为稠密矩阵

对于稀疏矩阵，函数cholinc( )来做不完全Cholesky分解
函数cholinc( )的另一个优点是它可用于计算实半正定矩阵

R=cholinc(X,DROPTOL)：
其中X和R的含义与函数chol( )中的变量相同，DROPTOL为不完全Cholesky分解的丢失容限。当DROPTOL设为0时，退化为完全Cholesky分解
R=cholinc(X,OPTS)：
其中OPTS为结构体，它有3个属性，即DROPTOL、MICHOL、RDIAG
1.DROPTOL为不完全Cholesky分解的丢失容限
2.MICHOL为1时采用改进算法的不完全Cholesky分解，否则不采用改进算法
3.RDIAG为1时R的对角元素中的零值替换成为DROPTPL的平方根，R=0时不做此替换
R=cholinc(X,‘0’)：完全Cholesky分解
[R,p]=cholinc(X,‘0’)：返回两个参数，并且不会返回出错信息。
当X是正定矩阵时，返回的上三角矩阵R满足$X=R^{\prime }\cdot R$，且p=0；
当X是非正定矩阵时，返回值p是正整数，R是上三角矩阵其阶数p-1，且满足X(1:p-1,1:p-1)=R’$\cdot$R
R=cholinc(X,‘inf’)：采用Cholesky-Infinity分解。此种分解是基于Cholesky分解的，但它可以处理实半正定矩阵

在MATLABR2018a中函数cholinc已经被删除，用ichol()函数替换

稀疏矩阵Cholesky分解示例：

1.2.1 一般方阵的高斯消去法分解（LU分解）

LU分解可将任意一个方阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即$A=LU$ (Lower Matrix，Upper Matrix)

[L,U]=lu(X)：
X是一个方阵,L是一个"心理"下三角矩阵,U为上三角矩阵,满足$X=L\cdot U$
[L,U,P]=lu(X)：
X是一个方阵,L是一个下三角矩阵,U为上三角矩阵,P为置换矩阵,满足$P\cdot X=L\cdot U$
Y=lu(X)：
X是一个方阵,把上三角矩阵和下三角矩阵合并在矩阵Y中，矩阵Y的对角元素为上三角矩阵的对角元素，即$Y=L+U-I$，置换矩阵P的信息丢失

考虑线性方程组$Ax=b$其中对矩阵A可以做LU分解，使得$A=L\cdot U$这样的线性方程组可以写成$L\cdot U\cdot x=b$，由于左除算符( \ )可以快速处理三角矩阵，因此可快速解出：$$x=U\backslash (L\backslash b)$$

利用LU分解来计算行列式的值和矩阵的逆

det(A)=det(L)*det(U)
inv(A)=inv(U)*inv(L)

进行LU分解示例：

对于稀疏矩阵，函数luinc( )来做不完全LU分解
(Incomplete LU decomposition)

[L,U]=luinc(X,DROPTPL)：
其中X、L、U的含义与函数lu()中的变量相同，DROPTPL为不完全LU分解的丢失容限。当DROPTPL设为0时，退化为完全LU分解
[L,U]=luinc(X,OPTS)：
其中OPTS为结构体有4个属性，即DROPTOL、MICHOL、RDIAG、THRESH
1.DROPTOL为不完全Cholesky分解的丢失容限
2.MICHOL为1时采用改进算法的不完全Cholesky分解，否则不采用改进算法
3.RDIAG为1时R的对角元素中的零值替换成为DROPTPL的平方根，R=0时不做此替换
4.THRESH是绕对角线旋转因子，其取值范围是[0,1]，当THRESH为0时强制绕对角线旋转，THRESH默认值是1
[L,U,P]=luinc(X,‘0’)：0级不完全LU分解
[L,U]=luinc(X,‘0’)：0级不完全LU分解
Y=luinc(X,‘0’)：0级完全LU分解

1.2.2 矩形矩阵的正交分解（QR分解）

QR分解把一个$m\times n$的矩阵A分解为一个正交矩阵Q和一个上三角矩阵R的乘积，即$$A=Q\cdot R$$

[Q,R]=qr(A)：其中矩阵R为与矩阵A具有相同大小的上三角矩阵，Q为正交矩阵，它们满足$A=Q\cdot R$。该调用方式适用于满矩阵和稀疏矩阵
[Q,R]=qr(A,0)：为"经济"方式的QR分解
设矩阵A是一个$m\times n$的矩阵
若$m＞n$，则只计算矩阵Q的前n列元素,R为$n\times n$矩阵
若$m\le n$，则与[Q,R]=qr(A)效果一致。该调用方式适用于满矩阵和稀疏矩阵
[Q,R,E]=qr(A)：R是上三角矩阵，Q为正交矩阵，E为置换矩阵
三者满足的关系：$A\cdot E=Q\cdot R$
程序选择一个合适的矩阵E使得abs(diag( R ))是降序排列的。该调用方式适用于满矩阵
[Q,R,E]=qr(A,0)："经济"方式的QR分解，其中E是一个置换矢量。
它们满足A(:,E)=$Q\cdot R$，该调用方式适用于满矩阵
R=qr(A)：返回上三角矩阵R，这里R=chol($A^{\prime }\cdot A$)。该调用方式适用于稀疏矩阵
R=qr(A,0)：以"经济"方式返回上三角矩阵R
[C,R]=qr(A,B)：
其中矩阵B必须与矩阵A具有相同的行数，矩阵R是上三角矩阵,$C=Q^{\prime }\cdot B$

通过QR分解分析矩阵的秩示例：

1.2.3 舒尔分解

舒尔分解定义式$$A=U\cdot S\cdot U^{\prime }$$其中A必须是一个方阵(Square Matrix)，
U是一个酉矩阵（或幺正矩阵)（图中共轭转置符号$†$ dagger,匕首）
S是一个块对角化矩阵,由对角线上的1×1和2×2块组成
特征值(eigenvalue)可以由矩阵S的对角块给出，而矩阵U给出比特征向量更多的数值特征。
此外，对缺陷矩阵也可以进行舒尔分解

[U,S]=schur(A)：返回酉矩阵U和块对角矩阵S
S=schur(A)：仅返回块对角矩阵
schur(A,‘real’)：返回的实特征值放在对角线上，而把复特征值放在对角线上的2×2块中
schur(A,‘complex’)：返回的矩阵S是上三角矩阵，并且如果矩阵A有复特征值，则矩阵S是复矩阵（指的是元素中含有复数的矩阵）
rsf2csf( )：把实数形式的舒尔矩阵转换成复数形式的舒尔矩阵

1.3 特征值和特征向量

1.3.0 特征值和特征向量的定义

假设A是一个$n\times n$的矩阵，A的特征值问题就是找到下面方程组的解：
$$A\cdot V=\lambda \cdot V$$其中，$\lambda$为标量，$V$为矢量

若把矩阵A的n个特征值放在矩阵P的对角线上，相应的特征向量按照特征值对应的顺序排列，作为矩阵V的列，特征值问题可改写为：$$A\cdot V=V\cdot D$$如果$V$是非奇异的，该问题可认为是一个特征值分解问题，关系式如下：$$A=V\cdot D\cdot V^{-1}$$

广义特征值问题是指方程$A\cdot X=\lambda \cdot B\cdot x$的非平凡解问题，其中A、B都是$n\times n$的矩阵，$\lambda$是标量。满足此方程的$\lambda$为广义特征值，对应的向量$x$为广义特征向量
如果$X$是一个列向量 $\vec{a}$ 的特征向量的矩阵，并且它的秩为n,那么特征向量线性相关。
如果不是这样，则称矩阵为缺陷阵
如果$X^{\prime }\cdot X=I$，则特征向量正交，这对于对称矩阵是成立的

1.3.1 特征值和特征向量的相关函数

eig(A)：求包含矩阵A的特征值的向量
[X,D]=eig(A)：产生一个矩阵A的特征值在对角线上的对角矩阵D和矩阵X，它们的列是相应的特征向量，满足$A\cdot X=X\cdot D$，为了得到有更好条件特征值的矩阵，要进行相似变换
[T,B]=balance(A)：找到一个相似变换矩阵T和矩阵B,使得它们满足$B=T-A\cdot T$，B是用balance命令求得的平衡矩阵
eig(A,‘nobalance’)：不经过平衡处理求得矩阵A的特征值和特征向量，也就是不进行平衡相似变换
eigs(A)：返回一个由矩阵A的部分特征值组成的向量，和eig命令一样，但不返回全部的特征值。如果不带有参量，则计算出最大的特征值，当计算所有特征值时，如果矩阵A的秩不小于6，则计算出6个特征值
eigs(f,n)：求出矩阵A的部分特征值。在使用一个矩阵列的线性运算符时，字符串f中包含的是M文件的文件名，n指定问题的阶次。用这种方法来求特征值比开始就用运算符来求要快
eigs(A,B,k,sigma)
求矩阵A的部分特征值,矩阵B的大小和A相同；
如果没有给出B=eye(size(A)),那么k就是要计算的特征值的个数；
如果k没有给出，就用小于6的数或者A的秩
sigma是一个实数或复数的移位参数，或下列文本字符串中的一个
文本字符串指明的是特征值的属性：
lm：为最大的特征值
sm：为最小的特征值
lr：为最大的实数部分
sr：最小的实数部分
be：为同时求的最大和最小的实数部分
condeig(A)：返回一个由矩阵A的特征值条件数组成的向量
[V,D,s]=condeig(A)：返回[V,D]=eig(A)和s=condeig(A)

1.3.2 特征值和特征向量的计算