YGC

12. Singular Value Decomposition (SVD)

3 min read

Write a Python function that approximates the Singular Value Decomposition on a 2x2 matrix by using the jacobian method and without using numpy svd function, i mean you could but you wouldn’t learn anything. return the result in this format.

Example:

  • Input:a = [[2, 1], [1, 2]]
  • Output:(array([[-0.70710678, -0.70710678],                       [-0.70710678,  0.70710678]]),       array([3., 1.]),       array([[-0.70710678, -0.70710678],              [-0.70710678,  0.70710678]]))

这道题是SVD（奇异值分解），一个矩阵分解成3个矩阵：

A=UΣVT

这里U和V都是正交矩阵。所以可以很容易地推出来：

AAT=UΣ2ATA=VΣ2

那么：

  • U就是AAT矩阵的特征向量所组成的矩阵，而Σ2就是相应的特征值。
  • V就是ATA矩阵的特征向量所组成的矩阵，而Σ2就是相应的特征值。

Σ2是一样的，求Σ就是对特征值进行开方。于是我们就可以求出来这三个矩阵。

分别去解U和V，到时候UΣVT一乘，它可能会不等于A，因为一个特征向量，它乘以一个数（可以是负数，改变方向），它还是一个特征向量。这就会出现问题。

为了解决这个问题，我就求V，然后再依赖V来求U，这样就能解决这个问题。

简单地推导一下，就可以得到U的计算公式，代码里求Σ的逆矩阵，我用了pinv，这样代码安全一点，因为不是所有矩阵都是inverse的，用pseudo inverse可以不出错。

import numpy as np   
def svd_2x2_singular_values(A: np.ndarray) -> tuple:  
    At = A.T  
    AtA = np.dot(At, A)  
    AAt = np.dot(A, At)  
    eig, V = np.linalg.eigh(AtA)  
    #eig, U = np.linalg.eigh(AAt)  
    S = np.sqrt(eig)  
  
    i = np.argsort(-S)  
    #U = U[:, i]  
    S = S[i]  
    V = V[:, i]  
    VT = V.T  
  
    S_inv = np.linalg.pinv(np.diag(S))  
    U = A.dot(V).dot(S_inv)  
  
    return (U, S, VT)

代码的中间对特征值进行排序，通常是排序的，然后我们可以扔掉比较小的特征值，这样可以减少运算，而且丢了可能还是好的，因为那些小的，可能是噪音。

如果不自己写代码，而直接调用svd()的话，那就是一行代码的事情：

np.linalg.svd(A)

Graph View