复现经典：《统计学习方法》第15章 奇异值分解

# 实现奇异值分解， 输入一个numpy矩阵，输出 U, sigma, V
# https://zhuanlan.zhihu.com/p/54693391

import numpy as np


#基于矩阵分解的结果，复原矩阵
def rebuildMatrix(U, sigma, V):
    a = np.dot(U, sigma)
    a = np.dot(a, np.transpose(V))
    return a


#基于特征值的大小，对特征值以及特征向量进行排序。倒序排列
def sortByEigenValue(Eigenvalues, EigenVectors):
    index = np.argsort(-1 * Eigenvalues)
    Eigenvalues = Eigenvalues[index]
    EigenVectors = EigenVectors[:, index]
    return Eigenvalues, EigenVectors


#对一个矩阵进行奇异值分解
def SVD(matrixA, NumOfLeft=None):
    #NumOfLeft是要保留的奇异值的个数，也就是中间那个方阵的宽度
    #首先求transpose(A)*A
    matrixAT_matrixA = np.dot(np.transpose(matrixA), matrixA)
    #然后求右奇异向量
    lambda_V, X_V = np.linalg.eig(matrixAT_matrixA)
    lambda_V, X_V = sortByEigenValue(lambda_V, X_V)
    #求奇异值
    sigmas = lambda_V
    sigmas = list(map(lambda x: np.sqrt(x)
                      if x > 0 else 0, sigmas))  #python里很小的数有时候是负数
    sigmas = np.array(sigmas)
    sigmasMatrix = np.diag(sigmas)
    if NumOfLeft == None:
        rankOfSigmasMatrix = len(list(filter(lambda x: x > 0,
                                             sigmas)))  #大于0的特征值的个数
    else:
        rankOfSigmasMatrix = NumOfLeft
    sigmasMatrix = sigmasMatrix[0:rankOfSigmasMatrix, :]  #特征值为0的奇异值就不要了

    #计算右奇异向量
    X_U = np.zeros(
        (matrixA.shape[0], rankOfSigmasMatrix))  #初始化一个右奇异向量矩阵，这里直接进行裁剪
    for i in range(rankOfSigmasMatrix):
        X_U[:, i] = np.transpose(np.dot(matrixA, X_V[:, i]) / sigmas[i])

    #对右奇异向量和奇异值矩阵进行裁剪
    X_V = X_V[:, 0:NumOfLeft]
    sigmasMatrix = sigmasMatrix[0:rankOfSigmasMatrix, 0:rankOfSigmasMatrix]
    #print(rebuildMatrix(X_U, sigmasMatrix, X_V))

    return X_U, sigmasMatrix, X_V

A = np.array([[1, 1, 1, 2, 2], [0, 0, 0, 3, 3], [0, 0, 0, 1, 1], [1, 1, 1, 0, 0],
              [2, 2, 2, 0, 0], [5, 5, 5, 0, 0], [1, 1, 1, 0, 0]])

A

array([[1, 1, 1, 2, 2],
       [0, 0, 0, 3, 3],
       [0, 0, 0, 1, 1],
       [1, 1, 1, 0, 0],
       [2, 2, 2, 0, 0],
       [5, 5, 5, 0, 0],
       [1, 1, 1, 0, 0]])

X_U, sigmasMatrix, X_V = SVD(A, NumOfLeft=3)

X_U

array([[ 1.96602638e-01, -5.12980706e-01, -6.20066911e-09],
       [ 3.08997616e-02, -8.04794293e-01,  1.69140901e-09],
       [ 1.02999205e-02, -2.68264764e-01,  5.63803005e-10],
       [ 1.76002797e-01,  2.35488225e-02, -7.63159275e-09],
       [ 3.52005594e-01,  4.70976451e-02, -1.52631855e-08],
       [ 8.80013984e-01,  1.17744113e-01, -3.81579637e-08],
       [ 1.76002797e-01,  2.35488225e-02, -7.63159275e-09]])

sigmasMatrix

array([[9.81586105e+00, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 5.25821946e+00, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00, 1.16381789e-07]])

X_V

array([[ 5.75872999e-01,  4.12749590e-02,  8.16496581e-01],
       [ 5.75872999e-01,  4.12749590e-02, -4.08248290e-01],
       [ 5.75872999e-01,  4.12749590e-02, -4.08248290e-01],
       [ 5.05512944e-02, -7.05297502e-01,  3.28082013e-17],
       [ 5.05512944e-02, -7.05297502e-01,  3.28082013e-17]])

# rebuild from U, sigma, V

rebuildMatrix(X_U, sigmasMatrix, X_V)

array([[ 1.00000000e+00,  1.00000000e+00,  1.00000000e+00,
         2.00000000e+00,  2.00000000e+00],
       [ 5.39915464e-17,  7.72260438e-16, -7.54662738e-16,
         3.00000000e+00,  3.00000000e+00],
       [ 9.57429619e-18,  2.48997260e-16, -2.59977132e-16,
         1.00000000e+00,  1.00000000e+00],
       [ 1.00000000e+00,  1.00000000e+00,  1.00000000e+00,
         1.25546281e-17,  1.25546281e-17],
       [ 2.00000000e+00,  2.00000000e+00,  2.00000000e+00,
         2.51092563e-17,  2.51092563e-17],
       [ 5.00000000e+00,  5.00000000e+00,  5.00000000e+00,
         9.74347659e-18,  9.74347659e-18],
       [ 1.00000000e+00,  1.00000000e+00,  1.00000000e+00,
         1.38777878e-17,  1.38777878e-17]])

from PIL import Image
import requests
from io import BytesIO

url = 'https://images.mulberry.com/i/mulberrygroup/RL5792_000N651_L/small-hampstead-deep-amber-small-classic-grain-ayers/small-hampstead-deep-amber-small-classic-grain-ayers?v=3&w=304'
response = requests.get(url)
img = Image.open(BytesIO(response.content))
img