「matlab学习笔记」MATLAB矩阵处理
中国大学MOOC 科学计算与MATLAB语言（点击此处跳转）
2.1 特殊矩阵
通用性的特殊矩阵
| 通用特殊矩阵函数 | 函数解释 |
|---|---|
| zeros函数 | 产生全0矩阵，即零矩阵 |
| ones函数 | 产生全1矩阵，即幺矩阵 |
| eye函数 | 产生对角线为1的矩阵。当矩阵是方阵时，得到一个单位矩阵 |
| rand函数 | 产生（0，1）区间均匀分布的随机矩阵 |
| randn函数 | 产生均值为0，方差为1的标准正态分布随机矩阵 |
以zeros函数为例，有3中调用格式：
- zeros(m)：产生m×m零矩阵。
- zeros(m,n)：产生m×n零矩阵。
- zeros(size(A))：产生与矩阵A同样大小的零矩阵。
>> A=zeros(2,3)
A =
0 0 0
0 0 0
>> zeros(size(reshape(A,3,2)))
ans =
0 0
0 0
0 0
拓展：
1.如何得到 \([a,b]\) 区间上均匀分布的随机整数
\(fix(a + (b-a+1) * rand(n))\)
\(\\\)
2.如何得到均值为 \(μ\)、方差为 \(σ^2\) 的随机数
\(μ + σ * randn(n)\)
%首先产生5阶两位随机整数矩阵A，再产生均值为0.6、方差为0.1的5阶正态分布随机矩阵B，最后验证(A+B)I=IA+BI（I为单位矩阵）
>> A=fix(10+(99-10+1)*rand(5));
>> B=0.6+sqrt(0.1)*randn(5);
>> C=eye(5);
>> (A+B)*C==C*A+B*C
ans =
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
用于专门学科的特殊矩阵
魔方矩阵（Magic Square）
- \(n\) 阶魔方阵由 \(1,2,3,…,n^2\) 共 \(n^2\) 个整数组成，且每行、每列以及主、副对角线上各 \(n\) 个元素之和都相等。
- \(n\) 阶魔方阵每行每列元素的和为 \(\cfrac {(1+2+3+…+ n^2)}{n}=\cfrac {(n+n^3)}{2}\)。
- \(n>2\) 时有很多不同的 \(n\) 阶魔方阵，MATLAB函数magic(n)产生一个特定的魔方阵。
>> M=magic(3)
M =
8 1 6
3 5 7
4 9 2
%产生8阶魔方阵，求其每行每列元素的和
>> M=magic(8);
>> sum(M(1,:))
ans =
260
>> sum(M(:,1))
ans =
260
范德蒙（Vandermonde）矩阵
对于向量 \(v=[v_1,v_2,…,v_n]\)，范得蒙矩阵的一般形式为：
- 在MATLAB中，函数vander(V)生成以向量V为基础的范得蒙矩阵。
- 范德蒙矩阵常用在各种通信系统的纠错编码中，例如，常用的Reed-Solomon编码即以范德蒙矩阵为基础。
>> A=vander(1:5)
A =
1 1 1 1 1
16 8 4 2 1
81 27 9 3 1
256 64 16 4 1
625 125 25 5 1
希尔伯特（Hilbert）矩阵
n阶希尔伯特矩阵的一般形式为：
- 希尔伯特矩阵的元素为 \(H(i,j)=1/(i+j-1)\)
- 在MATLAB中，生成n阶希尔伯特矩阵的函数是hilb(n)。
- 希尔伯特矩阵是著名的病态矩阵，即任何一个元素发生较小的变动，整个矩阵的值和逆矩阵都会发生很大变化。病态程度和矩阵的阶数相关，随着阶数的增加病态越严重。
>> format rat
>> H=hilb(4)
H =
1 1/2 1/3 1/4
1/2 1/3 1/4 1/5
1/3 1/4 1/5 1/6
1/4 1/5 1/6 1/7
伴随矩阵
设多项式 \(p(x)\) 为 \(a_nx^n+a_{n-1}x^{n-1}+…+a_1x+a_0\)，则多项式的伴随矩阵是：
- \(p(x)\) 称为 \(A\) 的特征多项式，方程 \(p(x)=0\) 的根称为 \(A\) 的特征值。
- MATLAB生成伴随矩阵的函数是compan(p)，其中 p 是一个多项式的系数向量，高次幂系数排在前，低次幂排在后。例如，生成多项式 \(x^3-2x^2-5x+6\) 的伴随矩阵。
- 可以求出伴随矩阵的特征值，该特征值等于多项式方程的根。
>> p=[1,-2,-5,6];
>> A=compan(p)
A =
2 5 -6
1 0 0
0 1 0
证明矩阵A的特征值等于多项式方程的根（大概率有误，如有问题欢迎指正）
按第n列展开行列式
记 \(A_n=|\lambda E-A|\)，则有
即
对比
又因为 \(a_n\neq0\)，所以 \(p(x)=0\) 即 \(A_n=|\lambda E-A|=0\)，得证
帕斯卡（Pascal）矩阵
- 根据二项式定理，\((x+y)^n\) 展开后的系数随着n的增大组成一个三角形表，这个三角形称为杨辉三角形。
- 把二项式系数依次填写在矩阵的左侧对角线上，然后提取左侧的n行n列元素即为n阶帕斯卡矩阵。
- 帕斯卡矩阵的第一行元素和第一列元素都为1，其余位置的元素是该元素的左边元素与上面元素相加，即 \(P(i,j)=P(i,j-1)+P(i-1,j)\)，且 \(P(i,1)=1,P(1,j)=1\)。
- 函数pascal(n)生成一个n阶帕斯卡矩阵。
%生成5阶帕斯卡矩阵，验证它的逆矩阵的所有元素也为整数
>> format rat
>> P=pascal(5)
P =
1 1 1 1 1
1 2 3 4 5
1 3 6 10 15
1 4 10 20 35
1 5 15 35 70
>> inv(P)
ans =
5 -10 10 -5 1
-10 30 -35 19 -4
10 -35 46 -27 6
-5 19 -27 17 -4
1 -4 6 -4 1
2.2 矩阵变换
对角阵
| 名词 | 解释 |
|---|---|
| 对角阵 | 只有对角线上有非零元素的矩阵 |
| 数量矩阵 | 对角线上的元素相等的对角矩阵 |
| 单位矩阵 | 对角线上的元素都为1的对角矩阵 |
（1）提取取矩阵的对角线元素
- diag(A)：提取矩阵A主对角线元素，产生一个列向量。
- diag(A,k)：提取矩阵A第k条对角线的元素，产生一个列向量。
（2）构造对角阵
- diag(V)：以向量 V为主对角线元素，产生对角矩阵。
- diag(V,k)：以向量 V为第k条对角线元素，产生对角矩阵。
1.对角阵左乘矩阵A，相当于将A的各行元素分别乘以对角阵的对角线元素。
2.对角阵右乘矩阵A，相当于将A的各列元素分别乘以对角阵的对角线元素。
记忆方法：将对角阵的位置看作二维数组的维度,在左边相当于行，在右边相当于列。
三角阵
| 名词 | 解释 |
|---|---|
| 上三角阵 | 矩阵的对角线以下的元素全为零的矩阵 |
| 下三角阵 | 对角线以上的元素全为零的矩阵 |
（1）上三角矩阵
- triu(A)：提取矩阵A的主对角线及以上的元素。
- triu(A,k)：提取矩阵A的第k条对角线及以上的元素。
>> triu(ones(4),-1)
ans =
1 1 1 1
1 1 1 1
0 1 1 1
0 0 1 1
（2）下三角矩阵
在MATLAB中，提取矩阵A的下三角矩阵的函数是tril，其用法与triu函数完全相同。
矩阵的转置
- 转置运算符是小数点后面接单引号（.'）。
- 共轭转置，其运算符是单引号（'），它在转置的基础上还要取每个数的复共轭。
如果矩阵的元素是实数，那么转置和共轭转置的结果是一样的。
>> A=[1,3;3+4i,1-2i]
A =
1.0000 + 0.0000i 3.0000 + 0.0000i
3.0000 + 4.0000i 1.0000 - 2.0000i
>> A.'
ans =
1.0000 + 0.0000i 3.0000 + 4.0000i
3.0000 + 0.0000i 1.0000 - 2.0000i
>> A'
ans =
1.0000 + 0.0000i 3.0000 - 4.0000i
3.0000 + 0.0000i 1.0000 + 2.0000i
矩阵的旋转
rot90(A,k)：将矩阵A逆时针方向旋转\(90^\circ\)的k倍，当k为1时可省略。
>> A=[1,3,2;-3,2,1;4,1,2]
A =
1 3 2
-3 2 1
4 1 2
>> rot90(A)
ans =
2 1 2
3 2 1
1 -3 4
>> rot90(A,2)
ans =
2 1 4
1 2 -3
2 3 1
矩阵的翻转
对矩阵实施左右翻转是将原矩阵的第一列和最后一列调换，第二列和倒数第二列调换，…，依此类推。
- fliplr(A)：对矩阵A实施左右翻转。
- flipud(A)：对矩阵A实施上下翻转。
注意区分矩阵的转置、共轭转置、旋转和翻转
%验证魔方阵的主对角线、副对角线元素之和相等
>> A=magic(5);
>> D1=diag(A);
>> sum(D1)
ans =
65
>> B=flipud(A);
>> D2=diag(B);
>> sum(D2)
ans =
65
矩阵求逆
- inv(A)：求方阵A的逆矩阵。
用求逆矩阵的方法解线性方程组
在线性方程组 \(Ax=b\) 两边各左乘 \(A^{-1}\)，得 \(x=A^{-1}b\)
>> A=[1,2,3;1,4,9;1,8,27];
>> b=[5;-2;6];
>> x=inv(A)*b
x =
23.0000
-14.5000
3.6667
>> x=A\b
x =
23.0000
-14.5000
3.6667
2.3 矩阵求值
矩阵的行列式值
det(A)：求方阵A所对应的行列式的值。
>> format rat
>> A=[1,3,2;-3,2,1;4,1,2]
A =
1 3 2
-3 2 1
4 1 2
>> det(inv(A))
ans =
1/11
>> 1/det(A)
ans =
1/11
矩阵的秩
- 矩阵线性无关的行数或列数称为矩阵的秩。
- rank(A)：求矩阵A的秩。
拓展：魔方阵的秩的规律
- 奇数阶魔方阵秩为 \(n\)，即奇数阶魔方阵是满秩矩阵
- 一重偶数阶魔方阵秩为 \(\cfrac {n}{2}+2\)（\(n\) 是 \(2\) 的倍数，但非 \(4\) 的倍数）
- 双重偶数阶魔方阵秩均为 \(3\)（阶数是 \(4\) 的倍数）
%求3~20阶魔方阵的秩
for n=3:20
r(n)=rank(magic(n));
end
bar(r)
grid on
axis([2,21,0,20])
[3:20;r(3:20)]
矩阵的迹
- 矩阵的迹等于矩阵的对角线元素之和，也等于矩阵的特征值之和。
- trace(A)：求矩阵A的迹。
>> A=[1,3,2;-3,2,1;4,1,2]
A =
1 3 2
-3 2 1
4 1 2
>> b = trace(A)
b =
5
>> t = sum(diag(A))
t =
5
矩阵的范数
矩阵或向量的范数用来度量矩阵或向量在某种意义下的长度。
（1）向量的范数
- 向量\(1\)——范数：向量元素的绝对值之和
- 向量\(2\)——范数：向量元素绝对值的平方和的平方根
- 向量\(\infty\)——范数：所有向量元素绝对值中的最大值
| 向量范数函数 | 函数解释 |
|---|---|
| norm(V)或norm(V,2) | 计算向量V的\(2\)——范数 |
| norm(V,1) | 计算向量V的\(1\)——范数 |
| norm(V,inf) | 计算向量V的\(\infty\)——范数 |
（2）矩阵的范数
- 矩阵A的\(1\)——范数：所有矩阵列元素绝对值之和的最大值
- 矩阵A的\(2\)——范数：A'A矩阵的最大特征值的平方根
- 矩阵A的\(\infty\)——范数：所有矩阵行元素绝对值之和的最大值
%MATLAB提供了求3种矩阵范数的函数
%其函数调用格式与求向量的范数的函数完全相同
>> x=[2 0 1;-1 1 0;-3 3 0]
x =
2 0 1
-1 1 0
-3 3 0
>> n = norm(x)
n =
4.7234
>> n = norm(x,1)
n =
6
矩阵的条件数
- 矩阵A的条件数等于A的范数与A的逆矩阵的范数的乘积。
- 条件数越接近于1，矩阵的性能越好，反之，矩阵的性能越差。
| 矩阵条件数函数 | 函数解释 |
|---|---|
| cond(A)或cond(A,2) | 计算A的\(2\)——范数下的条件数 |
| cond(A,1) | 计算A的\(1\)——范数下的条件数 |
| cond(A,inf) | 计算A的\(\infty\)——范数下的条件数 |
%求2~10阶希尔伯特矩阵的条件数
for n=2:10
c(n)=cond(hilb(n));
end
format long
c'
%运行结果
1.0e+13 *
0
0.000000000001928
0.000000000052406
0.000000001551374
0.000000047660725
0.000001495105864
0.000047536735631
0.001525757554777
0.049315394619572
1.602490962516758
%随着阶数的增加，希尔伯特矩阵的条件数不断增大，矩阵性能变差。
2.4 矩阵的特征值与特征向量
求矩阵的特征值与特征向量
在MATLAB中，计算矩阵的特征值和特征向量的函数是eig，常用的调用格式有两种：
- E=eig(A)：求矩阵A的全部特征值，构成向量E。
- [X,D]=eig(A)：求矩阵A的全部特征值，构成对角阵D，并产生矩阵X，X各列是相应的特征向量。
例1 设
又设 \(\lambda_i\) 为 \(R\) 的特征值，\(\lambda_j\) 为 \(S\) 的特征值，\(x_i=(\alpha_1,\alpha_2,\alpha_3)'\) 是 \(R\) 对应于 \(\lambda_i\) 的特征向量，\(y_j=(\beta_1,\beta_2)'\) 是 \(S\) 对应于 \(\lambda_j\) 的特征向量，试验证：
（1）\(\lambda_i\)、\(\lambda_j\)为 \(A\) 的特征值。
（2）\(p_i=(\alpha_1,\alpha_2,\alpha_3,0,0)'\) 是 \(A\) 对应于 \(\lambda_i\) 的特征向量，\(p_j=(0,0,0,\beta_1,\beta_2)'\) 是 \(A\) 对应于 \(\lambda_j\) 的特征向量。
R=[-1,2,0;2,-4,1;1,1,-6];
S=[1,2;2,3];
A=[R,zeros(3,2);zeros(2,3),S];
[X1,d1]=eig(R)
[X2,d2]=eig(S)
[X3,d3]=eig(A)
拓展：特征值的几何意义
设 \(A=\begin{bmatrix}3.8&0.6\\0.6&2.2\end{bmatrix}\)，其特征向量有 \(x_1=\begin{bmatrix}3\\1\end{bmatrix}\)，\(x_2=\begin{bmatrix}-1\\3\end{bmatrix}\)，对应的特征值分别为 \(\lambda_1=4\) 和 \(\lambda_2=2\)，令 \(y_1=Ax_1=\lambda_1x_1\)，\(y_2=Ax_2=\lambda_2x_2\)，我们讨论 \(y_1\) 与 \(x_1\) ，\(y_2\) 与 \(x_2\) 之间的关系。
更进一步，连续取单位向量x，让它大小保持为1，那么Ax就将四分之一圆弧进行拉伸，变成四分之一椭圆。
MATLAB提供了一个eigshow命令，可以演示向量x和Ax之间的关系。用鼠标拖动绿色的单位向量x绕原点转动，图中同步出现蓝色的Ax向量。Ax的大小在变化，方向也在变化，而且Ax的方向与x不一定相同。在变化过程中，x与Ax共线的位置称为特征方向。在特征方向上有Ax等于λx。
例2 已知大写字母M的各个结点坐标如表所示（第一行代表横坐标，第二行代表纵坐标）。
（1）绘制 M 的图形。
（2）设 \(A=\begin{bmatrix}1&0.5\\0&1\end{bmatrix}\)，用 A 对 M 的结点坐标进行变换，并绘制变换后的图形。
x=[0,0.5,0.5,3,5.5,5.5,6,6,3,0;0,0,6,0,6,0,0,8,1,8];
A=[1,0.5;0,1];
y=A*x;
subplot(2,2,1); %选择1号子图，详见专题四
fill(x(1,:),x(2,:),'r'); %绘制M的图形，并用红色（red）填充
subplot(2,2,2); %选择2号子图
fill(y(1,:),y(2,:),'r'); %绘制变换后的M图形，并用红色填充
定义变换矩阵A，再利用A对x进行变换，得到y矩阵，最后分别绘制变换前后的图形，M原来是正体，变换后改为斜体。
启示：在构建字库时，不必单独创建斜体字库，而只需对正体字库进行适当的线性变换即可，这样可以大大节省存储空间。
2.5 稀疏矩阵
矩阵的存储方式
- 完全存储方式：将矩阵的全部元素按列存储。
- 稀疏存储方式：只存储矩阵的非零元素的值及其位置，即行号和列号。
注意：采用稀疏存储方式时，矩阵元素的存储顺序并没有改变，也是按列的顺序进行存储。
稀疏存储方式的产生
（1）完全存储方式与稀疏存储方式之间的转化
- A=sparse(S)：将矩阵S转化为稀疏存储方式的矩阵A。
- S=full(A)：将矩阵A转化为完全存储方式的矩阵S。
（2）直接建立稀疏存储矩阵
- sparse(m,n)：生成一个m×n的所有元素都是零的稀疏矩阵。
- sparse(u,v,S)：其中u、v、S是3个等长的向量。S是要建立的稀疏存储矩阵的非零元素，u(i)、v(i)分别是S(i)的行和列下标。
>> A=sparse([1,2,2],[2,1,4],[4,5,-7])
A =
(2,1) 5
(1,2) 4
(2,4) -7
>> B=full(A)
B =
0 4 0 0
5 0 0 -7
（3）使用spconvert函数直接建立稀疏存储矩阵
B=spconvert(A)
其中，A 为一个 \(m\times3\) 或 \(m\times4\) 的矩阵，其每行表示一个非零元素，\(m\) 是非零元素的个数。
- A(i,1)表示第i个非零元素所在的行
- A(i,2)表示 第i个非零元素所在的列
- A(i,3)表示第i个非零元素值的实部
- A(i,4)表示第i个非零元素值的虚部
若矩阵的全部元素都是实数，则无须第4列
>> A=[2,2,1;2,1,-1;2,4,3]
A =
2 2 1
2 1 -1
2 4 3
>> B=spconvert(A)
B =
(2,1) -1
(2,2) 1
(2,4) 3
（4）单位矩阵的稀疏存储
speye(m,n)返回一个m×n的稀疏存储单位矩阵。
>> speye(3)
ans =
(1,1) 1
(2,2) 1
(3,3) 1
带状稀疏矩阵的稀疏存储
- 稀疏矩阵有两种基本类型：无规则结构的稀疏矩阵与有规则结构的稀疏矩阵。
- 带状稀疏矩阵就是一种十分典型的具有规则结构的稀疏矩阵，它是指所有非零元素集中在对角线上的矩阵。
| 带状稀疏矩阵函数 | 函数解释 |
|---|---|
| [B,d]=spdiags(A) | 从带状稀疏矩阵A中提取全部非零对角线元素赋给矩阵B及其这些非零对角线的位置向量d |
| A=spdiags(B,d,m,n) | 产生带状稀疏矩阵的稀疏存储矩阵A，其中m、n为原带状稀疏矩阵的行数与列数，矩阵B的第i列即为原带状稀疏矩阵的第i条非零对角线，向量d为原带状稀疏矩阵所有非零对角线的位置 |
>> A =[11,0,0,12,0,0;0,21,0,0,22,0;0,0,31,0,0,32;41,0,0,42,0,0;0,51,0,0,52,0]
A =
11 0 0 12 0 0
0 21 0 0 22 0
0 0 31 0 0 32
41 0 0 42 0 0
0 51 0 0 52 0
>> [B,d]=spdiags(A)
B =
0 11 12
0 21 22
0 31 32
41 42 0
51 52 0
d =
-3
0
3
%利用带状稀疏矩阵非零对角线元素组成的矩阵B，以及对角线位置组成的向量d，命令执行后产生一个稀疏存储矩阵A
>> A=spdiags(B,d,5,6)
A =
(1,1) 11
(4,1) 41
(2,2) 21
(5,2) 51
(3,3) 31
(1,4) 12
(4,4) 42
(2,5) 22
(5,5) 52
(3,6) 32
带状稀疏矩阵小结
利用 A=spdiags(B,d,m,n) 产生带状稀疏矩阵的稀疏存储
其中，m、n 为原带状矩阵的行数与列数，B 为 r×p 矩阵，这里 r=min(m,n)，p 为原带状矩阵所有非零对角线的条数，矩阵 B 的第 i 列即为原带状矩阵的第 i 条非零对角线。
取值方法：
- 若非零对角线上元素个数等于 r，则取全部元素；
- 若非零对角线上元素个数小于 r，则应该用零补足到 r 个元素。
补零原则：
- 若行数 < 列数，则主对角线以下在前面补 0，主对角线以上在后面补 0；
- 当行数 ≥ 列数，则主对角线以下在后面补 0，主对角线以上在前面补 0。
稀疏矩阵的应用实例
求下列三对角线性方程组的解
>> kf1=[1;1;2;1;0];
>> k0=[2;4;6;6;1];
>> k1=[0;3;1;4;2];
>> B=[kf1,k0,k1];
>> d=[-1;0;1];
>> A=spdiags(B,d,5,5);
>> f=[0;3;2;1;5];
>> x=A\f
x =
-0.1667
0.1111
2.7222
-3.6111
8.6111
注意区分稀疏矩阵和稀疏存储矩阵
友链
详细版笔记请跳转爆豆的博客
