Transcript Slides
MATLAB 程式設計:進階篇
線性代數
張智星 (Roger Jang)
[email protected]
http://mirlab.org/jang
台大資工系 多媒體檢索實驗室
MATLAB 程式設計進階篇:線性代數
反矩陣
反矩陣:
一個矩陣 A 的反矩陣可表示成
A1,滿足下列恆等式:
AA1 I
A1 A I
1
只有在 A 為方陣時, A 才存在。
1
若 A 不存在,則 A 稱為 Singular
MATLAB 程式設計進階篇:線性代數
反矩陣 (2)
inv:
MATLAB 的 inv 指令可用於計算反矩陣
範例6-1:inv01.m
A = pascal(4);
% 產生 4x4 的 Pascal 方陣
B = inv(A)
結果:
B=
4.0000 -6.0000
4.0000 -1.0000
-6.0000 14.0000 -11.0000
3.0000
4.0000 -11.0000 10.0000 -3.0000
-1.0000
3.0000 -3.0000 1.0000
Note that det(pascal(n)) is always equal to 1.
MATLAB 程式設計進階篇:線性代數
反矩陣 (3)
inv:
若矩陣 A 為 Singular (即其反矩陣不存在),則在使
用 inv 指令時,會產生警告訊息
範例6-2:inv02.m
A = [1 2 3; 4 5 6; 7 8 9];
B = inv(A)
結果:
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.541976e-018.
> In inv02 at 2
B=
1.0e+016 *
-0.4504 0.9007 -0.4504
0.9007 -1.8014 0.9007
-0.4504 0.9007 -0.4504
MATLAB 程式設計進階篇:線性代數
行列式
det:
欲計算矩陣 A 的行列式,可用 det 指令
範例6-3:det01.m
A = [1 3 4; -3 -4 -1; 2 2 5];
d = det(A)
結果:
d=
29
MATLAB 程式設計進階篇:線性代數
反矩陣公式
由 Crammer Rule 可知矩陣 A 的行列式和反矩陣有下列
關係式:
adj ( A)
A
A * adj ( A) A * I
A
1
其中 A 代表 A 的行列式,adj ( A) 代表 A 的 Adjoint
Matrix
Equivalent statements
A is singular
A1 does not exist
A 0
MATLAB 程式設計進階篇:線性代數
Cofactor and Adjoint
Take a 3x3 matrix A for example:
a11 a12
A a21 a22
a31 a32
a13
a23
a33
A11
cofactor( A) A21
A31
A12
A22
A32
adj A cofactor( A)T
A * adj ( A) A * I
a22
a
A13 32
a12
A23
a
A33 32
a12
a22
a23
a33
a13
a33
a13
a23
a21 a23
a31 a33
a11 a13
a31 a33
a
a
11 13
a21 a23
A a11 A11 a12 A12 a13 A13 ...
a21 a22
a31 a32
a11 a12
a31 a32
a11 a12
a21 a22
MATLAB 程式設計進階篇:線性代數
How to Prove the Formula?
How to prove A * adj ( A) A * I ?
a11 a12 a13 A11 A21 A31
A * adj ( A) a21 a22 a23 A12 A22 A32
a31 a32 a33 A13 A23 A33
a11 A11 a12 A12 a13 A13 a11 A21 a12 A22 a13 A23
a21 A11 a22 A12 a23 A13 a21 A21 a22 A22 a23 A23
a31 A11 a32 A12 a33 A13 a31 A21 a32 A22 a33 A23
A
0
0
0
A
0
0
0 A *I
A
a11 A31 a12 A32 a13 A33
a21 A31 a22 A32 a23 A33
a31 A31 a32 A32 a33 A33
Note that :
A a11 A11 a12 A12 a13 A13
a21 A21 a22 A22 a23 A23
a31 A31 a32 A32 a33 A33
MATLAB 程式設計進階篇:線性代數
Quiz Candidates
If A is a 3x3 matrix, what is adj(A)?
a11 a12
A a21 a22
a31 a32
a13
a23 adj A ?
a33
And can you demonstrate that
adj A* A A * I
MATLAB 程式設計進階篇:線性代數
反矩陣與行列式之驗證 (1/2)
若 A 為整數矩陣,則 A 乘上 A1 必為整數矩陣,可
驗証如下:
範例6-4:det02.m
A = [1 3 4; -3 -4 -1; 2 2 5];
det(A)*inv(A)
結果:
ans =
-18
13
2
-7
-3
4
13
-11
5
MATLAB 程式設計進階篇:線性代數
反矩陣與行列式之驗證 (2/2)
若 A 為整數矩陣,將 inv(A) 以有理形式(Rational
Format,即分子和分母都是整數的分數)來表示,即
可察覺出它和行列式的關係
從這裡可以很明顯的看出,inv(A) 中每個元素的分母
值,就是 det(A)
範例6-5:det03.m
A = [1 3 4; -3 -4 -1; 2 2 5];
format rat
% 以有理形式表示數值
inv(A)
format short
% 改回預設的數值表示形式
ans =
結果:
-18/29
-7/29
13/29
13/29
-3/29
-11/29
2/29
4/29
5/29
MATLAB 程式設計進階篇:線性代數
固有值與固有向量
一個方陣 A 的固有向量(Eigenvector) x 與
固有值 (Eigenvalue) 的關係式如下:
Ax x
上式可化簡成
A I ) x 0
由於 x 不是一個零向量,因此 A I 必須是
Singular,上式才會有解。當 A I是 Singular
時,其行列式為零:
A I 0
MATLAB 程式設計進階篇:線性代數
固有值分解 (1/2)
若 A 為 n×n 的矩陣,則上式為 的 n 次多項式 ,代表 將
有 n 個解 1 n ,滿足下列關係式:
Ax1 1 x1
Ax x
n n
n
或可寫成矩陣形式: AX XD
其中
| |
X x1 xn
| |
1 0 0
D 0 0
0 0 n
如果 X 1 存在,則由上式可得矩陣A的固有值分解
(Eigenvalue decomposition):
A XDX 1
MATLAB 程式設計進階篇:線性代數
固有值分解 (2/2)
If A is symmetric, then X
Suppose A is 3x3, then
1
XT
A XDX 1 XDX T
|
x1
|
|
x2
|
| 1
x3 0
| 0
0
2
0
0
0
3
| 1 x1T
T
x2 x3 2 x2
|
| 3 x3T
1 x1 x1T 2 x2 x2T 3 x3 x3T
|
x1
|
|
x1T
x2T
x3T
MATLAB 程式設計進階篇:線性代數
計算固有值
eig:
MATLAB 的 eig 指令可用於計算矩陣的固有值及固有
向量。若只想計算固有值時,可輸入如下:
範例6-6:eig01.m
A = magic(5);
lambda = eig(A)
結果:
lambda =
65.0000
-21.2768
-13.1263
21.2768
13.1263
MATLAB 程式設計進階篇:線性代數
計算固有值與固有向量
若要同時計算固有值及固有向量,須提供兩個輸出引數
範例6-7:eig02.m
A = magic(5);
[X, D] = eig(A)
結果:
X=
-0.4472
-0.4472
-0.4472
-0.4472
-0.4472
0.0976
0.3525
0.5501
-0.3223
-0.6780
-0.6330
0.6780
-0.2619
0.5895
0.3223
-0.1732
-0.3915
-0.5501
0.3915
0.1732
-0.3525
-0.5895
0.2619
-0.0976
0.6330
D=
65.0000
0
0
0 -21.2768
0
0
0 -13.1263
0
0
0
0
0
0
0
0
0
21.2768
0
0
0
0
0
13.1263
請驗證看看
X*D*inv(X)是否等於A。
MATLAB 程式設計進階篇:線性代數
固有向量和固有值的展示
Try “eigshow” to plot the trajectory of a
linear transform in 2D
MATLAB 程式設計進階篇:線性代數
奇異值與奇異向量
一個矩陣 A 與其奇異值(Singular Value) 及奇異向
量(Singular Vectors)u 與 v 之間存在下列的關係式:
Av u
t
A u v
若將所有的行向量 u 並排成矩陣 U,所有的行向量 v 並排成矩
陣 V。則上式可寫成:
AV U
T
A U V
其中 Σ 的主對角線即是對應的 值,其餘元素為零
MATLAB 程式設計進階篇:線性代數
奇異值分解
若 A 的維度是 m×n,則 U、Σ、V 的維度分別是 m×m、m×n
以及 n×n。一般而言,U 和 V 均是 Unitary 矩陣(即矩陣內
的行向量均兩兩相互垂直,且行向量的長度均為 1),滿足
下列條件:
UU T I
T
VV I
因此矩陣 A 可寫成:
A U V
T
上式稱為奇異值分解(Singular Value Decomposition)
MATLAB 程式設計進階篇:線性代數
計算奇異值與奇異向量
svd:
svd 指令可用於計算矩陣的奇異值及奇異向量
範例6-8:svd01.m
A = [3 5; 4 7; 2 1; 0 3];
請驗證看看
[U, S, V] = svd(A)
U*S*V’是否等於A。
結果:
U=
S=
V=
-0.5577
0.0881
-0.6954
0.4447
10.4517
0
-0.7714
0.0333
0.2489
-0.5848
0
-0.1771
0.6471
0.5453
0.5025
0
0
-0.2504
-0.7565
0.3965
0.4558
0
0
1.9397
-0.4892
0.8722
-0.8722
-0.4892
MATLAB 程式設計進階篇:線性代數
計算奇異值與奇異向量 (2)
若 A 為 m×n 的矩陣且 m >> n,則我們可在原先的
svd 指令加入另一個輸入引數 0,使其產生的 U 及 S
矩陣具有較小的維度
範例6-9:svd02.m
A = [3 5; 4 7; 2 1; 0 3];
[U, S, V] = svd(A, 0)
結果:
U=
-0.5577
-0.7714
-0.1771
-0.2504
0.0881
0.0333
0.6471
-0.7565
S=
10.4517
0
V=
-0.4892
-0.8722
0
1.9397
0.8722
-0.4892
請驗證看看
U*S*V’是否等於A。
MATLAB 程式設計進階篇:線性代數
線性聯立方程式
聯立線性方程式:
線性代數中最重要的問題,即是解決聯立線性方程式。
一組線性方程式可用矩陣表示如下:
AX = B
其中 A、B 是已知矩陣,而 X 則是未知矩陣
MATLAB 程式設計進階篇:線性代數
線性聯立方程式之可能情況
線性聯立方程式:
為簡化起見,我們可以假設 A、X、B 的維度分別是
m×n、n×1、m×1,其中 m 代表方程式的數目,n 則是
未知數的數目,可分成三種情況:
若 m = n,則方程式的個數和未知數的個數相等,通常會有一
組解滿足 AX = B
若 m > n,則方程式的個數大於未知數的個數,通常無一解可
滿足 AX = B,但我們可轉而求取最小平方解(Least-Squares
2
Solution)X,使得 AX B 為最小值
若 m < n,則方程式的個數小於未知數的個數,通常有無限多
組解可滿足 AX = B,我們可尋求一基本解(Basic Solution)X,
使得 X 最少包含 m-n 個零元素
MATLAB 程式設計進階篇:線性代數
斜線和反斜線運算
斜線和反斜線運算:
MATLAB 提供一個反斜線運算(Back Slash Operator,
即「\」)使得 X=A\B 能滿足上述三種情況
反斜線運算又稱「左除」(Left Division)
MATLAB 也提供了斜線運算(Slash Operator,即「/」)
斜線運算又稱「右除」(Right Division)
可用於對付 XA = B 的方程組。
MATLAB 程式設計進階篇:線性代數
斜線和反斜線運算:記憶法
整理:MATLAB 常用的數學函數
方程式形式
MATLAB
解法
AX = B
X = A\B(左除)
XA = B
X = B/A(右除)
在上表中,欲解 AX = B 或 XA = B,我們可以想像在
等號兩邊各除以 A,並依 A 的位置分別取用「左除」或
「右除」
MATLAB 程式設計進階篇:線性代數
範例:通過三點的二次曲線
斜線和反斜線運算:
範例6-10:leftDiv01.m
A = vander(1:3);
B = [6; 11; 18];
X = A\B
結果:
X=
1.0000
2.0000
3.0000
>> A*X – B
Ans =
0
0
0
此例代表通過 (1, 6), (2, 11), (3, 18) 的二次曲線為
y x2 2x 3
MATLAB 程式設計進階篇:線性代數
範例:最小平方解
斜線和反斜線運算:
當 m > n 時,「左除」可以找到最小平方解
範例6-11:leftDiv02.m
A = [2 -1; 1 -2; 1 1];
B = [2; -2; 1];
X = A\B
結果:
X=
1.0000
1.0000
在上例中,我們有 3 個方程式,但卻只有 2 個未知數,
此 3 個方程式在 x-y 平面並未交於一點,故嚴格地說,
此方程組無解,而 MATLAB「左除」找到的 X 為最小平
2
方解,使得 AX B 為最小
MATLAB 程式設計進階篇:線性代數
範例:基本解
斜線和反斜線運算:
當 m < n 時,「左除」可以找到基本解
範例6-12:leftDiv03.m
A = [1 2 3; 4 5 6];
B = [7; 8];
X = A\B
結果:
X=
-3.0000
0
3.3333
基本解至少有 n-m 個零。
MATLAB 程式設計進階篇:線性代數
使用「左除」進行最小平方法
問題:在平面上找出一點 P,使得 P 到下列三條
直線的距離平方和為最小:
3x - 4y = 2
4x y = - 3
x + y =1
Hint : The shortest distance between
a point x0 , y0 and a line ax by c 0
is
解法:
2
ax0 by0 c
a b
2
2
2
.
2
3x - 4y - 2 4x y 3 x y 1
L ( x, y )
A
x
b
, where
5
17
2
4/5
3/ 5
2/5
x
A 4 / 17 1 / 17 , b 3 / 17 , x
y
1 / 2 1 / 2
1 / 2
2
MATLAB 程式設計進階篇:線性代數
類似問題
解析幾何
問題:在平面上找出一
點 P,使得 P 到下列X的
Y為最小
X:三條直線、三點
Y:距離和、距離平方和
古典幾何
問題:在平面上找出一
點 P,使得 P 到三角型
的X的Y為最小
X:三邊、三頂點
Y:距離和、距離平方和
MATLAB 程式設計進階篇:線性代數
本章指令彙整
與線性代數相關的函式,彙整如下:
類別
函式
功能
矩陣
相關性質
det
norm
normest
null
行列式
矩陣或向量的 norm
估測矩陣的 2-norm
Null 空間
orth
垂直基底
(Orthonormal Basis)
rank
矩陣的 rank
rref
Reduced Row Echelon
Form
subspace
子空間的夾角
trace
矩 陣 對角線元素的總
和
MATLAB 程式設計進階篇:線性代數
本章指令彙整 (2)
與線性代數相關的函式,彙整如下:
類別
函式
矩陣函式
expm
funm
logm
sqrtm
功能
矩陣的指數函數
矩陣的一般函數
矩陣的對數函數
矩陣的開平方
左除及右除,用於解線
「\」及「/」
性方程式
Cholesky 分解
chol
不完全的 Cholesky 分
choline
解
MATLAB 程式設計進階篇:線性代數
本章指令彙整 (3)
與線性代數相關的函式,彙整如下:
類別
函式
矩陣函式
cond
condest
Inv
lu
luinc
qr
lscov
nnls
功能
矩 陣 的 Condition
Number,以測試其接
近 Singular 的程度
1-norm
Condition
Number 的估計
反矩陣
LU 分解
不完全 LU 分解
QR 分解
QR 分解
非負的最小平方法