1 滿秩分解（Full-Rank Factorization）

• 滿秩分解定理

  設$m\times n$矩陣$A$的秩為$r>0$，則存在$m\times r$矩陣$B$（列滿秩矩陣）和$r\times n$矩陣$C$（行滿秩矩陣）使得
  $$A=BC$$
  并且$rank(B)=rank(C)=r$

  滿秩分解不唯一

  定理：設$A$為$m\times n$矩陣，且$rank(A)=r$，存在$m$階可逆矩陣$P$和$n$階可逆矩陣$Q$，使得$A=P\left[\begin{array}{cc}{I}_r& 0\\ 0& 0\end{array}\right]Q$。

  證明滿秩分解定理：

  $\because \left[\begin{array}{cc}{I}_r& 0\\ 0& 0\end{array}\right]=\left[\begin{array}{c}{I}_r\\ 0\end{array}\right]\left[\begin{array}{cc}{I}_r& 0\end{array}\right]$

  $\therefore A=P\left[\begin{array}{c}{I}_r\\ 0\end{array}\right]\left[\begin{array}{cc}{I}_r& 0\end{array}\right]Q$

  則令$P\left[\begin{array}{c}{I}_r\\ 0\end{array}\right]=B$，$\left[\begin{array}{cc}{I}_r& 0\end{array}\right]Q=C$，可得到$A=BC$

  $\because$ $P,C$是可逆矩陣，$B$的$r$個列是$P$的前$r$列；$C$的$r$個行是$Q$的前$r$行

  $\therefore$ $rank(B)=rank(C)=r$

• 滿秩分解步驟

  1. 設$A$為$m\times n$矩陣，首先求$rank(A)$
  2. 取$A$的$j_1,j_2,...j_r$列構成$B_{m\times r}$
  3. 取$A$的Hermite標準型（即行最簡行矩陣）$H$的前$r$行構成矩陣$C_{r\times n}$
  4. 則$A=BC$就是矩陣$A$的一個滿秩分解

• Python求解滿秩分解

  import numpy as np
  from sympy import Matrix'''
  Full-Rank Factorization
  @params: A Matrix
  @return: F, G Matrix
  '''
  def full_rank(A):r = A.rank()A_arr1 = np.array(A.tolist())# 求解A的最簡行階梯矩陣，要轉換成list，再轉換成arrayA_rref = np.array(A.rref()[0].tolist())k = [] # 存儲被選中的列向量的下標# 遍歷A_rref的行for i in range(A_rref.shape[0]):# 遍歷A_rref的列for j in range(A_rref.shape[1]):# 遇到1就說明找到了A矩陣的列向量的下標# 這些下標的列向量組成F矩陣，然后再找下一行if A_rref[i][j] == 1:k.append(j)break# 通過選中的列下標，構建F矩陣       B = Matrix(A_arr1[:,k])# G就是取行最簡行矩陣A的前r行構成的矩陣C = Matrix(A_rref[:r])return B, Cif __name__ == "__main__":# 表示矩陣AA = np.array([[1, 1, 0], [0, 1, 1], [-1, 0, 0], [1, 1, 1]])A = Matrix(A)B, C = full_rank(A)print("B:", B)print("C:", C)
  

2 三角分解（Triangular Factorization）

• $LU$分解定義

  如果有一個矩陣$A$，我們能表示下三角矩陣$L$和上三角矩陣$U$的乘積，稱為$A$的三角分解或$LU$分解。

  

  更進一步，我們希望下三角矩陣的對角元素都為$1$

  

• $LU$分解定理

  若$A$是==$n$階非奇異矩陣==，則存在唯一的單位下三角矩陣$L$和上三角矩陣$U$使得$A=LU$的充分必要條件是$A$的所有順序主子式均非零（這一條件保證了對角線元素非零），即${\Delta }_k\ne 0(k=1,...,n?1)$

• $LU$分解步驟

  設$A$為$n\times n$矩陣

  1. 進行初等行變換（注意：不涉及行交換的初等變換），從第$1$行開始，到第$n$行結束。將第$i$行第$i$列以下的元素全部消為$0$
  2. 這樣操作后得到的矩陣即為$U$
  3. 構造對角線元素全為$1$的單位下三角矩陣$L$，$L$的剩余元素通過構建方程組的形式來求解。

• Python求解$LU$分解

• $LU$分解的實際意義

  • 解線性方程組

    假設我們有一個線性方程組$Ax=b$，其中$A$是一個非奇異矩陣，而$b$是一個列向量。通過$LU$分解，我們可以將方程組轉化為兩個簡化的方程組$Ly=b$和$Ux=y$，其中$L$是下三角矩陣，$U$是上三角矩陣。這兩個方程組分別易于求解。

    具體：

    首先，通過前代法（forward substitution）解$Ly=b$，然后通過回代法（backward substitution）解$Ux=y$。這樣，我們就得到了方程組的解。

• $LDU$分解定理

  設$A$是==$n$階非奇異矩陣==，則存在唯一的單位下三角矩陣$L$，對角矩陣$D=diag(d_1,d_2,...,d_n)$和上三角矩陣$U$使得$A=LDU$的充分必要條件是$A$的所有順序主子式均非零（這一條件保證了對角線元素非零），即${\Delta }_k\ne 0(k=1,...,n?1)$并且$d_1=a_{11},d_k=\frac{{\Delta }_k}{{\Delta }_{k+1}},k=2,...,n$

• $LDU$分解步驟

  設$A$為$n\times n$矩陣

  1. 先求$LU$分解
  2. 將$U$的對角線元素提出來構成對角矩陣$D$
  3. $U$中的元素$u_{ij}$除以$d_i$，其中$d_i$表示第$i$個對角元素。這樣操作得到變換后的$U$

• Python求解$LDU$分解

  import numpy as np
  from sympy import Matrix
  import pprintEPSILON = 1e-8def is_zero(x):return abs(x) < EPSILONdef LU(A):# 斷言A必須是非奇異方陣Aassert A.rows == A.cols, "Matrix A must be a square matrix"assert A.det() != 0, "Matrix A must be a nonsingular matrix"n = A.rowsU = A# 構建出U矩陣# 將U轉換成list，再轉換成arrayU = np.array(U.tolist())# 遍歷U的每一行利用高斯消元法for i in range(n):# 判斷U[i][i]是否為0assert not is_zero(U[i][i]), "主元為0，無法進行LU分解"# 對i+1行到n行進行消元for j in range(i + 1, n):# 計算消元因子factor = U[j][i] / U[i][i]# 對第j行進行消元for k in range(i, n):U[j][k] -= factor * U[i][k]# 消元后的矩陣U則是最終U矩陣U = Matrix(U)# 根據LU = A，得到L矩陣L = A * U.inv()return L, Udef LDU(A):L, U = LU(A)D = Matrix(np.diag(np.diag(U)))U = D.inv() * Ureturn L, D, Uif __name__ == '__main__':A = np.array([[2, 3, 4], [1, 1, 9], [1, 2, -6]])A = Matrix(A)'''# test LU分解L, U = LU(A)pprint.pprint("L:")pprint.pprint(L)pprint.pprint("U:")pprint.pprint(U)'''# test LDU分解L, D, U = LDU(A)pprint.pprint("L:")pprint.pprint(L)pprint.pprint("D:")pprint.pprint(D)pprint.pprint("U:")pprint.pprint(U)
  

• $PLU$分解

  PLU 分解是將矩陣$A$分解成一個置換矩陣$P$、單位下三角矩陣$L$和上三角矩陣$U$的乘積，即
  $$A=PLU$$
  之前$LU$分解中限制了行交換，如果不可避免的必須進行行互換，我們就需要進行$PLU$分解。

  實際上只需要把$A=LU$變成$P^{?1}A=P^{?1}PLU$就可以了，實際上所有的$A=LU$都可以寫成$P^{?1}A=LU$的形式，由于左乘置換矩陣$P^{?1}$是在交換行的順序，所以由$P^{?1}A=P^{?1}PLU$推得適當的交換$A$的行的順序，即可將$A$ 做 $LU$ 分解。當$A$沒有行互換時，$P$就是單位矩陣。

  事實上，所有的方陣都可以寫成 $PLU$ 分解的形式，事實上，$PLU$ 分解有很高的數值穩定性，因此實用上是很好用的工具。

  有時為了計算上的方便，會同時間換行與列的順序，此時會將 $A$ 分解成
  $$A=PLUQ$$
  其中 $P$、$L$、$U$ 同上，$Q$ 是一個置換矩陣。

3 正交三角分解（QR Factorization）

• $QR$分解定理

  設$A$是$m\times n$實（復）矩陣，$m\ge n$且其$n$個列向量線性無關，則存在$m$階正交（酉）矩陣$Q$和$n階$非奇異實（復）上三角矩陣$R$使得
  $$A=Q\left[\begin{array}{c}R\\ 0\end{array}\right]$$

• $QR$分解步驟

  設$A$為$3\times 3$矩陣，即$A=({\alpha }_1,{\alpha }_2,{\alpha }_3)$。則：

  1. 正交化：${\beta }_1={\alpha }_1$，${\beta }_2={\alpha }_2?k_{21}{\beta }_1$，${\beta }_3={\alpha }_3?k_{31}{\beta }_1?k_{32}{\beta }_2$，其中$k_{21}=\frac{<{\alpha }_2,{\beta }_1>}{<{\beta }_1,{\beta }_1>}$，$k_{31}=\frac{<{\alpha }_3,{\beta }_1>}{<{\beta }_1,{\beta }_1>}$，$k_{31}=\frac{<{\alpha }_3,{\beta }_2>}{<{\beta }_2,{\beta }_2>}$。

  2. 單位化得到矩陣$Q$：$Q=(\frac{{\beta }_1}{||{\beta }_1||},\frac{{\beta }_2}{||{\beta }_2||},\frac{{\beta }_3}{||{\beta }_3||})$

  3. 計算得到矩陣$R$：
     $$\left(\begin{array}{ccc}||{\beta }_1||& & \\ & ||{\beta }_2||& \\ & & ||{\beta }_3||\end{array}\right)\left(\begin{array}{ccc}1& k_{21}& k_{31}\\ & 1& k_{32}\\ & & 1\end{array}\right)=\left(\begin{array}{ccc}||{\beta }_1||& ||{\beta }_1||\times k_{21}& ||{\beta }_1||\times k_{31}\\ & ||{\beta }_2||& ||{\beta }_2||\times k_{32}\\ & & ||{\beta }_3||\end{array}\right)$$

  4. 這樣，$A=QR$

• Python求解$QR$分解

  常規計算：

  import numpy as np
  import sympy
  from sympy import Matrix
  from sympy import *
  import pprint#正交三角分解（QR）
  a = [[1, 1, -1],[-1, 1, 1],[1, 1, -1],[1, 1, 1]]# a = [[1,1,-1],
  #                   [1,0,0],
  #                   [0,1,0],
  #                   [0,0,1]]
  A_mat = Matrix(a)#α向量組成的矩陣A
  # A_gs= GramSchmidt(A_mat)
  A_arr = np.array(A_mat)
  L = []
  for i in range(A_mat.shape[1]):L.append(A_mat[:,i])
  #求Q
  A_gs = GramSchmidt(L)#α的施密特正交化得到β
  A_gs_norm = GramSchmidt(L,True)#β的單位化得到vA = []for i in range(A_mat.shape[1]):for j in range(A_mat.shape[0]):A.append(A_gs_norm[i][j])#把數排成一行A_arr = np.array(A)
  A_arr = A_arr.reshape((A_mat.shape[0],A_mat.shape[1]),order = 'F')#用reshape重新排列（‘F’為豎著寫）
  #得到最后的Q
  Q = Matrix(A_arr)#求RC = []
  for i in range(A_mat.shape[1]):for j in range(A_mat.shape[1]):if i > j:C.append(0)elif i == j:t = np.array(A_gs[i])m = np.dot(t.T,t)C.append(sympy.sqrt(m[0][0]))else:t = np.array(A_mat[:,j])x = np.array(A_gs_norm[i])n = np.dot(t.T,x)
  #             print(n)C.append(n[0][0])
  # C_final為R          
  C_arr = np.array(C)
  # print(C_arr)
  C_arr = C_arr.reshape((A_mat.shape[1],A_mat.shape[1]))
  R = Matrix(C_arr)pprint.pprint("Q:")
  pprint.pprint(Q)
  pprint.pprint("R:")
  pprint.pprint(R)
  

  調用庫函數

  # 求矩陣A的QR分解，保留根號
  Q_, R_ = A_mat.QRdecomposition()
  pprint.pprint("Q_:")
  pprint.pprint(Q_)
  pprint.pprint("R_:")
  pprint.pprint(R_)
  assert Q_ == Q, "Q_ != Q"
  assert R_ == R, "R_ != R"
  

4 奇異值分解（SVD）

• $SVD$定理

  設$A$是$m\times n$矩陣，且$rank(A)=r$，則存在$m$階酉矩陣$U$和$n$階酉矩陣$V$使得
  $$A=U\left[\begin{array}{cc}\Sigma & 0\\ 0& 0\end{array}\right]V^H$$
  其中$\Sigma =diag({\sigma }_1,...,{\sigma }_r)$，且${\sigma }_1\ge ...\ge {\sigma }_r>0$。

  $\sigma$為$A$的奇異值，具體含義這里不在敘述，但需要記住的是${\sigma }^2$是$A^{H}A$的特征值，也是$A{A}^H$的特征值，且：

  1. $A^{H}A$與$A{A}^H$的特征值均為非負數
  2. $A^{H}A$與$A{A}^H$的非零特征值相同，并且非零特征值的個數（重特征值按重數計算）等于$rank(A)$

  所以我們求$\Sigma$就轉換成求這兩個矩陣其中一個的特征值。

• $SVD$分解步驟

  1. 求$A^{H}A$的$n$個特征值，即計算$|\lambda I?A^{H}A|=0$。得到特征值：${\lambda }_1,...,{\lambda }_r,{\lambda }_{r+1}=0,...,{\lambda }_n=0$，其中$r=rank(A)$。

  2. 將$r$個奇異值（即非零特征值開根號）從大到小排列組成對角矩陣，再添加額外的$0$構成${\Sigma }_{m\times n}$矩陣。
     $${\Sigma }_{m\times n}=\left(\begin{array}{ccccc}\sqrt[]{{\lambda }_1}& ...& 0& 0& 0\\ 0& \sqrt[]{{\lambda }_2}& 0& 0& 0\\ ...& ...& ...& ...& ...\\ ...& ...& ...& \sqrt[]{{\lambda }_r}& ...\\ 0& ...& 0& 0& 0\end{array}\right)$$

  3. 求特征值：${\lambda }_1,...,{\lambda }_r,{\lambda }_{r+1}=0,...,{\lambda }_n=0$對應的特征向量${\xi }_1,...,{\xi }_n$：當$\lambda ={\lambda }_1$時，$(\lambda I?A^{H}A)\times {\xi }_1=0$，解得${\xi }_1$，同理，計算其余特征向量。

  4. 因為${\xi }_1,...,{\xi }_n$相互正交，我們還需要進行單位化，得到$v_1,...,v_n$，即$v_1=\frac{{\xi }_1}{||{\xi }_1||},...,v_n=\frac{{\xi }_n}{||{\xi }_n||}$。則$V=(v_1,...,v_n)$。

  5. 根據$A=U_{m\times m}{\Sigma }_{m\times n}{V}_{n\times n}^H$，可得$U_1=A{V}_{n\times n}{\Sigma }_{r\times n}^{?1}$（注意，$\sigma$此時為${\Sigma }_{m\times n}$的前$r$行），易知$U_1$為$m\times r$的矩陣，我們還需要擴充$U_2$，其為$m\times (m?r)$矩陣。

  6. 取$U_1^H{U}_2=0$，取$U_2$，必須要單位化$U_2$，這樣$U=[U_1:U_2]$

  7. 則$A=U_{m\times m}{\left[\begin{array}{cc}\Sigma & 0\\ 0& 0\end{array}\right]}_{m\times n}{V}_{n\times n}^H$

• Python求解奇異值分解

  import numpy as np
  from sympy import Matrix
  import pprintA = np.array([[1,0],[0,1],[1,1]])
  # 求A的奇異值分解
  U, sigma, VT = np.linalg.svd(A)
  print ("U:", U)
  print ("sigma:", sigma)
  print ("VT:", VT)