matlab/Matlab数值线性代数

矩阵分析

函数概览

函数调用方式 解释
det(A) 计算A的行列式
diag(A) 取A对角元构成向量或利用向量A生成对角阵
expm(A) 计算A的矩阵指数运算
inv(A) 计算A的逆矩阵
rank(A) 计算A的秩
tril(A) 将A化为下三角矩阵(其余元素变为0)
[V D]=eig(A) 计算A的特征值保存到D,特征向量保存到V
[V J]=jordan(A) 计算A的jordan分解,$AV=VJ$
[U S V]=svd(A) 计算A的奇异值分解,$A=USV’$
[L U]=lu(A) 计算A的LU分解,$A=LU$
[Q R]=qr(A) 计算A的QR分解,$A=QR$

矩阵的范数

  1. norm(A):用于计算矩阵的$2-$范数(也称为谱范数),定义为:

  2. norm(A, 'fro'):用于计算矩阵的Frobenius范数,即把矩阵拉成一个长向量,再计算其2-向量范数,常用于矩阵的误差分析,定义为:

  3. norm(A, 1):用于计算矩阵的$1-$范数,也称列范数

    其中向量$x$的1-范数为向量$x$所有元素绝对值的和

  4. norm(A, inf):用于计算矩阵的$\infty-$范数,也称行范数

    其中向量$x$的$\infty$范数即$x$所有元素的最大绝对值

  5. sum(sum(A~=0.0)):计算矩阵的$0-$范数(注意这并不是真正的范数),含义为矩阵A的非零元素的个数,记为

    $0-$范数是对矩阵稀疏性的一种最直接、最准确的度量

  6. [U S V] = svd(A); norm(diag(S), 1):计算矩阵的核范数(矩阵A的所有奇异值之和),对核范数的约束是低
    秩逼近与低秩求解最常用的方法,定义为:

低秩约束实例

  1. 强制降秩去噪模型:通过获取截取主成分来降秩,例如

  2. 低秩惩罚项模型:通过对参数$λ>0$的大小调整,以及矩阵本身奇异值与0的逼近程度来自动控制降秩的幅度

    设$A^*$的奇异值为$\bold{\sigma}=[\sigma_1, \sigma_2, … ,\sigma_n]^T$,该模型还可以表示为:

  3. 松弛化低秩惩罚项模型:为了使模型具有更好的收敛性,可以将奇异值的$ℓ_0$范数更改为$ℓ_1$范数,对应以下模型

  4. 半正定核范数正则化问题,在 $𝑨∗$半正定时,其核范数等于迹:

矩阵标量特征的验证

同阶方阵或“内维”相等的矩阵相乘,交换次序迹不变
1
2
3
4
A = rand(3, 3); B = rand(3, 3); C = rand(3, 4); D = rand(4, 3);
tAB = trace(A*B); tBA = trace(B*A);
tCD = trace(C*D); tDC = trace(D*C);
disp([tAB, tBA, tCD, tDC])
经典代数结论$|AB| = |BA| = |A|·|B|$
1
2
3
4
d_A_B = det(A)*det(B);
d_AB = det(A*B); d_BA = det(B*A);
d_CD = det(C*D); d_DC = det(D*C);
disp([d_A_B, d_AB, d_BA, d_CD, d_DC])

需要注意到:内维相等的矩阵相乘交换次序行列式可能会改变

矩阵变换与特征值分解

  • [R,ci] = rref(A)实现初等变换将$A$化为行阶梯型矩阵,其中$R$储存了行阶梯型的矩阵,$ci$表示阶梯元素(线性独立)的列指标

  • null(A) 得到 A 的零空间(核)的一组标准正交基

  • orth(A)得到 A 的值空间(像)的一组标准正交基

rref 化矩阵为行阶梯型矩阵

1
2
A = magic(4);
[R, ci] = rref(A)
1
2
3
4
5
6
7
8
9
10
R =

1 0 0 1
0 1 0 3
0 0 1 -3
0 0 0 0

ci =

1 2 3

上述结果意为A的秩等于3,行阶梯型共有3层,阶梯元素在前三列

矩阵零空间及其含义
1
2
3
4
5
6
A = reshape(1:15,5,3);
X = null(A); % 零空间向量
S = A * X; % 代入后应该得到零向量
n=size(A,2); % A的列数
l=size(X,2); % X的列数(即零空间维数)
n-l==rank(A) % 判断矩阵秩是否等于(列数-零空间维数)

返回结果为logical 1

实矩阵的特征值分解
1
2
3
4
5
6
7
8
9
A=[1, -3 ;2, 2/3];
[V,D]=eig(A); % 该矩阵有两个复特征根(D对角元),V为列特征向量构成的矩阵
[VR,DR]=cdf2rdf(V,D); % 将复数对角型转换为实数块对角型(VR*DR/VR=A)
A1 = real(V*D/V);
A2 = VR*DR/VR;
disp([A1, A2]);
err1 = norm(A-A1, 'fro');
err2 = norm(A-A2, 'fro');
disp([err1, err2]);

Matlab解线性方程组

A\b是Matlab最高效的解线性方程组的內建函数,A\B可解矩阵方程。

Matlab对于矩阵A列满秩且有解的情况将提供唯一解,对A列不满秩且有解的情况将随机提供一个解。无解问题将提供最小二乘解。

1
2
3
4
5
6
7
8
9
A = reshape(1:12, 4, 3);
b = (13:16)';
ra = rank(A); rab = rank([A, b]);
disp([ra, rab]);
xs = A\b; %求出方程组一个特解记为xs
xg = null(A); %求对应齐次方程零空间(基础解系)
c = rand(1);
ba = A*(xs+c*xg); %检验任意特解+对应齐次特解是否符合方程组
norm(ba-b) %向量2-范数的误差分析,新版本误差变小
函数Inv与左右除
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
rng default
format
A = gallery('randsvd', 300, 2e13, 2);
x = ones(300, 1);
b = A * x;
cond(A);

tic
xi = inv(A)*b; % 先求逆,再相乘浪费时间
ti = toc; % 将当前计时存储在变量ti中


eri = norm(x-xi); % 计算解的绝对误差
rei = norm(A*xi-b)/norm(b); % 计算回代相对误差

disp([cond(A)]); % 2.0040e+13
disp([ti, eri, rei]) % 0.0016 0.1032 0.0048

线性方程组的数值迭代法

从最初的线性方程组Ax=b出发,首先令$A=M-N$,然后有$Mx=Nx+b$,即$x=M^{-1}Nx+M^{-1}b = Bx+f$

数值迭代法会将上述等式右边的x代入上一步的解向量,从而得出下一步的解向量,如:

这种算法在$||B||<1$的情况可保持收敛性,根据$M,N$的不同选取,可以衍生出不同的基本迭代算法。

雅可比迭代法

采用矩阵分裂记号:$A = D-L-U$,并有$M=D, N = L+U$,其中$D$为对角阵,$L$为下三角阵,$U$为上三角阵。雅可比迭代法的矩阵表现形式为

雅可比迭代法数值实现:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
rng default
A=rand(400,400)+5e2*eye(400);
%雅可比迭代法需要对角元明显大于其他元素,否则一旦||B||>=1就会使算法不收敛
x=rand(400,1);
b=A*x;
D=diag(diag(A)); %取出A的对角元
N = D-A; %N=L+U=D-A
B = D\N; f = D\b; %x=Bx+f雏形已确定

x_old = zeros(400,1);%迭代初值设为0向量
err = 1e12; %初始误差假设很大
iter = 0; %记录总迭代次数
while(err>1e-6) %当相邻两步绝对误差不超过1e-6迭代即结束
iter = iter + 1;
x_new = B*x_old + f;
err = norm(x_new-x_old);%记录该步结果与上一步结果的最新绝对误差
x_old = x_new;
end
disp(iter); %共迭代19次
err = norm(x_new-x) %最终绝对误差

变分法解线性方程组

定义$\phi(x) = \frac12x^TAx-b^Tx$,则线性方程组$Ax=b$的解集等于 $\text{argmin}_{x\in \mathbb{R}^n}\phi(x)$的解集,将线性方程组求解问题转化为最优化问题的方法,称为变分法

由于矩阵A的正定性,最优化问题$\text{argmin}_{x\in \mathbb{R}^n}\phi(x)$为一个向量空间的凸问题,凸问题仅有一个全局最优解,可以通过梯度下降法求解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
rng default,A0 = rand(300,300);
A=A0+A0'+1e3*eye(300); % A为一个随机的正定矩阵即可
% 在很多实际问题中,比雅可比,高斯-塞德尔方法更加可行
x=rand(300,1);
b=A*x;

x_old = zeros(300,1); % 初始迭代值仍设为0向量
err = 1e12;
iter = 0;

while(err>1e-6) % 相邻两步误差不超过1e-6则终止迭代
iter = iter + 1;
p = b - A*x_old; % 迭代方向(目标函数负梯度)
lambda = (p'*p)/((A*p)'*p); % 最速迭代步长
x_new = x_old + lambda*p;
err = norm(x_new-x_old);% 统计迭代当前步的误差
x_old = x_new;
end
err = norm(x_new-x) % 最终绝对误差 9.2710e-08
iter % 循环次数 8
  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!

扫一扫,分享到微信

微信分享二维码
  • Copyrights © 2020-2021 chenk
  • 由 帅气的CK本尊 强力驱动
  • 访问人数: | 浏览次数:

请我喝杯咖啡吧~

支付宝
微信