首页 模拟笔记

# 奇异值分解

# 特征值和特征向量

特征值和特征向量的定义如下

𝐴𝑥=𝜆𝑥

其中 A 是一个𝑛×𝑛的实对称矩阵，𝑥是一个𝑛维向量，则我们说𝜆是矩阵 A 的一个特征值，而𝑥是矩阵 A 的特征值𝜆所对应的特征向量。

# SVD 定义

D 也是对矩阵进行分解，但是和特征分解不同，SVD 并不要求要分解的矩阵为方阵。假设我们的矩阵 A 是一个𝑚×𝑛的矩阵，那么我们定义矩阵 A 的 SVD 为：

A=UΣVTA = U\Sigma V^T

其中 U 是一个𝑚×𝑚的矩阵，Σ 是一个𝑚×𝑛的矩阵，除了主对角线上的元素以外全为 0，主对角线上的每个元素都称为奇异值，V 是一个𝑛×𝑛的矩阵。U 和 V 都是酉矩阵，即满足UTU=I,VTV=IU^TU = I,V^TV = I

Untitled

如果我们将 A 的转置和 A 做矩阵乘法，那么会得到𝑛×𝑛的一个方阵ATAA^TA。既然ATAA^TA 是方阵，那么我们就可以进行特征分解，得到的特征值和特征向量

(𝐴^𝑇𝐴)𝑣_𝑖=𝜆_𝑖𝑣_𝑖

这样我们就可以得到矩阵ATAA^TA 的 n 个特征值和对应的 n 个特征向量 vv 了。将ATAA^TA 的所有特征向量张成一个n×nn \times n 的矩阵 V，那么AATAA^T 就是的所有特征向量就是矩阵UU，

Untitled

也就说 sig 矩阵可以用𝐴𝑣_𝑖=𝜎_𝑖u_i 𝜎_𝑖就是奇异值

Untitled

其中 SVD 的物理含义就是

Untitled

这就表明任意的矩阵 A 是可以分解成三个矩阵相乘的形式。V 表示了原始域的标准正交基，U 表示经过 A 变换后的 co-domain 的标准正交基，Σ 表示了 V 中的向量与 U 中相对应向量之间的关系。我们仔细观察上图发现，线性变换 A 可以分解为旋转、缩放、旋转这三种基本线性变换。

sig 是对角阵，表示奇异值，A 矩阵的作用是将一个向量在 V 这组正交基向量的空间旋转，并对每个方向进行了一定的缩放，缩放因子就是各个奇异值。然后在 U 这组正交基向量的空间再次旋转。可以说奇异值分解将一个矩阵原本混合在一起的三种作用效果，分解出来了。

# SVD 应用

接下来我们从分解的角度重新理解前面的表达式，我们把原来的矩阵 A 表达成了 n 个矩阵的和：

Untitled

其中奇异值就是每个矩阵的权重，越往前的比例越大。

我们用一个图像压缩的例子来说明。我们知道，电脑上的图像（特指位图）都是由像素点组成的，所以存储一张 1000×622 大小的图片，实际上就是存储一个 1000×622 的矩阵，共 622000 个元素。这个矩阵用 SVD 可以分解为 622 个矩阵之和，如果我们选取其中的前 100 个之和作为对图像数据的近似，那么只需要存储 100 个奇异值，100 个向量和 100 个向量，共计 100×(1+1000+622)=162300 个元素，大约只有原始的 26% 大小。

# 解 SVD

其中解 SVD 有几种办法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
struct Functions
{
	// -----------------------------
	
		
	float3 randomUnitVector()
	{
		float3 unnormalized;
		for (int i = 0; i < 3; ++i)
		{
			unnormalized[i] = RandBBSfloat(i);
		}

			
		return unnormalized;
	}

	float3 svd_1d(float3x3 A)
	{
		float epsilon = 1e-10;
		

		float3 x = randomUnitVector();
		float3 lastV = float3(0, 0, 0);
		float3 currentV = x;
			
		float3x3 B = mul(A, transpose(A));
			
		for(int i = 0;i<9;i++)
		{
			lastV = currentV;
			currentV = mul(B, lastV);
			currentV = normalize(currentV);

			if (abs(dot(currentV, lastV)) > 1 - epsilon)
			{
				// converged
				break;
			}
		}
		return currentV;
		
	}

	void svd(float3x3 A, out float3x3 us,out float3 singularValues,  out float3x3 vs )
	{
		
	
		float3 u;
		float3 v;
	    
		for (int i = 0; i < 3; ++i)
		{
			float3x3 matrixFor1D = A;

			for (int j = 0; j < i; ++j)
			{
					
				u = us[j];
				v = vs[j];
				float3x3 OuterProduct = float3x3(u.x * v.x, u.x * v.y, u.x * v.z,
					u.y * v.x, u.y * v.y, u.y * v.z,
					u.z * v.x, u.z * v.y, u.z * v.z);
					
				matrixFor1D -= singularValues[j] * OuterProduct;
			}

			float3 u = svd_1d(matrixFor1D);
			float3 v_unnormalized = mul(transpose(A), u);
			float sigma = length(v_unnormalized);
			float3 v = v_unnormalized / sigma;
				

			singularValues[i] = sigma;
			us[i] = u;
			vs[i] = v;
		}
	}
	void ConverMatrix(float3x3 Inmat, out float4x4 OutMat)
	{
		OutMat =  float4x4(Inmat[0].x,Inmat[0].y,Inmat[0].z,0,
			Inmat[1].x,Inmat[1].y,Inmat[1].z,0,
			Inmat[2].x,Inmat[2].y,Inmat[2].z,0,
			0,0,0,1);
		
	}
	void ConverMatrix(float4x4 Inmat,out float3x3 OutMat)
	{
		OutMat =  float3x3(Inmat[0].x,Inmat[0].y,Inmat[0].z,
			Inmat[1].x,Inmat[1].y,Inmat[1].z,
			Inmat[2].x,Inmat[2].y,Inmat[2].z);
	}
	
};

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
import time
import numpy as np
from numpy.linalg import norm

from random import normalvariate
from math import sqrt

def randomUnitVector(n):
    unnormalized = [normalvariate(0, 1) for _ in range(n)]
    theNorm = sqrt(sum(x * x for x in unnormalized))
    return [x / theNorm for x in unnormalized]

def svd_1d(A, epsilon=1e-10):
    ''' The one-dimensional SVD '''

    n, m = A.shape
    x = randomUnitVector(min(n,m))
    lastV = None
    currentV = x

    if n > m:
        B = np.dot(A.T, A)
    else:
        B = np.dot(A, A.T)

    iterations = 0
    while True:
        iterations += 1
        lastV = currentV
        currentV = np.dot(B, lastV)
        currentV = currentV / norm(currentV)

        if abs(np.dot(currentV, lastV)) > 1 - epsilon:
            print("converged in {} iterations!".format(iterations))
            return currentV

def svd(A, k=None, epsilon=1e-10):
    '''
        Compute the singular value decomposition of a matrix A
        using the power method. A is the input matrix, and k
        is the number of singular values you wish to compute.
        If k is None, this computes the full-rank decomposition.
    '''
    A = np.array(A, dtype=float)
    n, m = A.shape
    svdSoFar = []
    if k is None:
        k = min(n, m)

    for i in range(k):
        matrixFor1D = A.copy()

        for singularValue, u, v in svdSoFar[:i]:

            matrixFor1D -= singularValue * np.outer(u, v)

        if n > m:
            v = svd_1d(matrixFor1D, epsilon=epsilon)  # next singular vector
            u_unnormalized = np.dot(A, v)
            sigma = norm(u_unnormalized)  # next singular value
            u = u_unnormalized / sigma
        else:
            u = svd_1d(matrixFor1D, epsilon=epsilon)  # next singular vector
            v_unnormalized = np.dot(A.T, u)
            sigma = norm(v_unnormalized)  # next singular value
            v = v_unnormalized / sigma

        svdSoFar.append((sigma, u, v))

    singularValues, us, vs = [np.array(x) for x in zip(*svdSoFar)]
    return singularValues, us.T, vs

def QR_HouseHolder(mat: np.array):
    cols = mat.shape[1]

    Q = np.eye(cols)
    R = np.copy(mat)

    for col in range(cols - 1):
        a = np.linalg.norm(R[col:, col])
        e = np.zeros((cols - col))
        e[0] = 1.0
        num = R[col:, col] - a * e
        den = np.linalg.norm(num)

        u = num / den
        H = np.eye(cols)
        H[col:, col:] = np.eye((cols - col)) - 2 * u.reshape(-1, 1).dot(u.reshape(1, -1))
        R = H.dot(R)

        Q = Q.dot(H)

    return Q, R

def QR_GivenRot(mat: np.array):
    rows, cols = mat.shape

    R = np.copy(mat)
    Q = np.eye(cols)

    for col in range(cols):
        for row in range(col + 1, rows):
            if abs(R[row, col]) < 1e-6:
                continue

            f = R[col, col]
            s = R[row, col]
            den = np.sqrt(f * f + s * s)
            c = f / den
            s = s / den

            T = np.eye(rows)
            T[col, col], T[row, row] = c, c
            T[row, col], T[col, row] = -s, s

            R = T.dot(R)

            Q = T.dot(Q)

    return Q.T, R

# ---------------------------------------
def QR_Basic(mat: np.array):
    eig = np.copy(mat)

    rows, cols = mat.shape
    eigV = np.eye(cols)
    for _ in range(50):
        q, r = QR_HouseHolder(eig)
        eig = r.dot(q)
        eigV = eigV.dot(q)

    return np.diag(eig), eigV

# ---------------------------------------
def QR_Heisenberg(mat: np.array):
    rows, cols = mat.shape
    Hsbg = np.copy(mat)

    eigv = np.eye(cols)
    for row in range(1, rows - 1):
        array = mat[row:, row - 1]

        # householder
        a = np.linalg.norm(array)

        e = np.zeros((1, rows - row))[0]
        e[0] = 1.0
        num = array - a * e
        u = num / np.linalg.norm(num)

        H = np.eye(rows - row) - 2 * u.reshape(-1, 1).dot(u.reshape(1, -1))

        # Hsbg
        T = np.eye(rows)
        T[row:, row:] = H

        Hsbg = T.T.dot(Hsbg).dot(T)
        eigv = eigv.dot(T.T)

    eigv = eigv.T
    # Given
    for _ in range(50):
        q, r = QR_GivenRot(Hsbg)
        Hsbg = r.dot(q)
        eigv = eigv.dot(q)
    return np.diag(Hsbg), eigv

# -----------------------------
def Jacobi_Basic(mat: np.array):
    rows, cols = mat.shape

    eig = np.copy(mat)
    eigV = np.eye(rows, cols)

    for _ in range(300):
        maxRow, maxCol = 0, 0
        maxValue = -1
        for row in range(rows):
            for col in range(cols):
                if row == col:
                    continue
                if abs(eig[row, col]) > maxValue:
                    maxValue = abs(eig[row, col])
                    maxRow = row
                    maxCol = col
        a = 0
        if abs(eig[maxRow, maxRow] - eig[maxCol, maxCol]) < 1e-5:
            a = 0.25 * np.pi
        else:
            a = 0.5 * np.arctan((2 * eig[maxRow, maxCol]) / (eig[maxRow, maxRow] - eig[maxCol, maxCol]))

        P = np.eye(rows)
        P[maxRow, maxRow] = np.cos(a)
        P[maxCol, maxCol] = np.cos(a)
        P[maxRow, maxCol] = -np.sin(a)
        P[maxCol, maxRow] = np.sin(a)
        eig = P.T.dot(eig).dot(P)
        eigV = eigV.dot(P)

    return np.diag(eig), eigV

# ----------------
def SortEig(eig, eigV):
    e = np.copy(eig)
    v = np.copy(eigV)
    print(e);
    indices = np.argsort(e)
    print(indices);
    e = np.sort(e)

    for i, j in enumerate(indices):
        v[:, i] = eigV[:, j]

    return e, v

def SVD_Solver(mat, Func):
    aTa = mat.T.dot(mat)
    aaT = mat.T.dot(mat)
    print(f"Ata:{ aTa}" )
    eig, eigv = Func(aTa)

    eig, eigv = SortEig(eig, eigv)

    rows, cols = mat.shape
    U = eigv
    S = np.zeros_like(mat)
    for i in range(rows):
        S[i, i] = eig[i]

    eig, eigv = Func(aaT)
    eig, eigv = SortEig(eig, eigv)
    V = eigv

    return U, np.sqrt(S), V

def SVD1(mat: np.array, Method: int):
    # Method: 0: QR_Basic, 1:QR_Heisenberg, 2:Jacobi_Basic
    U, S, V = SVD_Solver(mat, Jacobi_Basic)

    return U, S, V

mat = np.array([[3.0, 2.0, 2.0],
                [2.0, 5.0, 1.0],
                [2.0, 1.0, 4.0], ])

start_time = time.perf_counter()

end_time = time.perf_counter()
usetime = end_time - start_time
# print(f"Method QR_Basic use: {usetime} " )
# print(u)
# print(s)
# print(v)
# print("\n")

start_time = time.perf_counter()
u, s, v = svd(mat)
end_time = time.perf_counter()
usetime = end_time - start_time
print(f"Method Jacobi_Basic use: {usetime}" )
print(s)
# print(s)
# print(v)
# print("\n")

start_time = time.perf_counter()
u, s, v = np.linalg.svd(mat)#SVD(mat, 0)
end_time = time.perf_counter()
usetime = end_time - start_time
print(f"Method Real use: {usetime}")

print(u)
print(s)
print(v)
# print("\n")

其中这段代码最快是 jacobi 算法