写了一下, 效率很低, 数据样本是从GSLA扒下来的, 有空继续改进

import numpy as np
import math

# Given a matrix and its transposed
matrix = np.array([[2, 2], [1, 1]])
matrix_transposed = np.matrix.transpose(matrix)

# Generate (A^T)A and A(A^T)
left_symmetric_matrix = np.matmul(matrix, matrix_transposed)
right_symmetric_matrix = np.matmul(matrix_transposed, matrix)

left_singular_values, left_singular_vectors = np.linalg.eig(
    left_symmetric_matrix)
right_singular_values, right_singular_vectors \
    = np.linalg.eig(right_symmetric_matrix)


# Sort the eigenvalues
U = np.array([ele for _, ele in sorted(
    zip(left_singular_values, left_singular_vectors), reverse=True)])

# Construct a diagonal
Sigma = np.zeros(matrix.shape)
np.fill_diagonal(Sigma, [math.sqrt(abs(val)) for val in left_singular_values])

# V^T
VT = np.array([ele for _, ele in
               sorted(zip(right_singular_values, right_singular_vectors), reverse=True)])
VT = np.matrix.transpose(VT)

print("U=", '\n', U, end='\n\n')
print("Sigma=", '\n', Sigma, end='\n\n')
print("V^T=", '\n', VT, end='\n\n')

@bintou

非常牛逼！很好。

@osmanthus 在 SVD 分解的py代码? 中说：

import numpy as np
import math

Given a matrix and its transposed

matrix = np.array([[2, 2], [1, 1]])
matrix_transposed = np.matrix.transpose(matrix)

Generate (A^T)A and A(A^T)

left_symmetric_matrix = np.matmul(matrix, matrix_transposed)
right_symmetric_matrix = np.matmul(matrix_transposed, matrix)

left_singular_values, left_singular_vectors = np.linalg.eig(
left_symmetric_matrix)
right_singular_values, right_singular_vectors
= np.linalg.eig(right_symmetric_matrix)

Sort the eigenvalues

U = np.array([ele for _, ele in sorted(
zip(left_singular_values, left_singular_vectors), reverse=True)])

Construct a diagonal

Sigma = np.zeros(matrix.shape)
np.fill_diagonal(Sigma, [math.sqrt(abs(val)) for val in left_singular_values])

V^T

VT = np.array([ele for _, ele in
sorted(zip(right_singular_values, right_singular_vectors), reverse=True)])
VT = np.matrix.transpose(VT)

print("U=", '\n', U, end='\n\n')
print("Sigma=", '\n', Sigma, end='\n\n')
print("V^T=", '\n', VT, end='\n\n')

@osmanthus

In [5]: run svd.py
File "/home/lbwang/SCode/PyCode/svd.py", line 31
print("U=", '\n', U, end='\n\n')
^
SyntaxError: invalid syntax

@bintou 要用py3跑

@osmanthus 折腾啊。。。Py3的numpy还没安装。慢慢来吧。。。

@bintou sudo pip3 install numpy

@osmanthus 搞定了。我问一下，你输出的是列向量还是行向量？

@bintou ?

@bintou 已经是矩阵的形式打印出来了, 每一列是一个列向量

import numpy as np
import math

# Given a matrix and its transposed
matrix = np.array([[2, 2], [1, 1]])
matrix_transposed = np.matrix.transpose(matrix)

# Generate (A^T)A and A(A^T)
left_symmetric_matrix = np.matmul(matrix, matrix_transposed)
right_symmetric_matrix = np.matmul(matrix_transposed, matrix)

left_singular_values, left_singular_vectors = np.linalg.eig(
    left_symmetric_matrix)
right_singular_values, right_singular_vectors \
    = np.linalg.eig(right_symmetric_matrix)


# Sort the eigenvalues
U = np.array([ele for _, ele in sorted(
    zip(left_singular_values, left_singular_vectors), reverse=True)])

# Construct a diagonal
Sigma = np.zeros(matrix.shape)
np.fill_diagonal(Sigma, [math.sqrt(abs(val)) for val in left_singular_values])

# V^T
VT = np.array([ele for _, ele in
               sorted(zip(right_singular_values, right_singular_vectors), reverse=True)])
VT = np.matrix.transpose(VT)

print("U=", '\n', U, end='\n\n')
print("Sigma=", '\n', Sigma, end='\n\n')
print("V^T=", '\n', VT, end='\n\n')

print("Back to A=", '\n', np.matmul(np.matmul(U, Sigma), VT))

这个加了一条将分解出来的矩阵乘回去的语句

@osmanthus
Sage有SVD的实现，比如：

m = matrix(RDF,[[2,2],[1,1]])
m.SVD()

得到的答案与你的有符号上的区别，不知道为什么。

@bintou 感觉应该分解不唯一吧, 乘回去都是一样的诶