Update README.md

README.md

@@ -171,3 +171,97 @@ print()
 print(Y[:, :, 1])
 print()
 ```
+
+<h2 align="center">低秩时序模型</h2>
+<p align="right"><a href="#从线性代数到张量分解"><sup>▴ 回到顶部</sup></a></p>
+
+```python
+import numpy as np
+
+def compute_rse(var, var_hat):
+    return np.linalg.norm(var - var_hat, 2) / np.linalg.norm(var, 2)
+
+def laplacian(T, tau):
+    ell = np.zeros(T)
+    ell[0] = 2 * tau
+    for k in range(tau):
+        ell[k + 1] = -1
+        ell[-k - 1] = -1
+    return ell
+
+def prox(z, w, lmbda, denominator):
+    T = z.shape[0]
+    temp = np.fft.fft2(lmbda * z - w) / denominator
+    temp1 = 1 - T / (lmbda * np.abs(temp))
+    temp1[temp1 <= 0] = 0
+    return np.fft.ifft2(temp * temp1).real
+
+def update_z(y_train, pos_train, x, w, lmbda, eta):
+    z = x + w / lmbda
+    z[pos_train] = (lmbda / (lmbda + eta) * z[pos_train] 
+                    + eta / (lmbda + eta) * y_train)
+    return z
+
+def update_w(x, z, w, lmbda):
+    return w + lmbda * (x - z)
+
+def laplacian_conv_2d(y_true, y, lmbda, gamma, eta, tau, maxiter = 50):
+    M, N = y.shape
+    pos_train = np.where(y != 0)
+    y_train = y[pos_train]
+    pos_test = np.where((y_true != 0) & (y == 0))
+    y_test = y_true[pos_test]
+    z = y.copy()
+    w = y.copy()
+    ell_1 = laplacian(M, tau)
+    ell_2 = laplacian(N, tau)
+    denominator = lmbda + gamma * np.fft.fft2(np.outer(ell_1, ell_2)) ** 2
+    del y_true, y
+    show_iter = 10
+    for it in range(maxiter):
+        x = prox(z, w, lmbda, denominator)
+        z = update_z(y_train, pos_train, x, w, lmbda, eta)
+        w = update_w(x, z, w, lmbda)
+        if (it + 1) % show_iter == 0:
+            print(it + 1)
+            print(compute_rse(y_test, x[pos_test]))
+            print()
+    return x
+```
+
+```python
+import numpy as np
+np.random.seed(1)
+from skimage import color
+from skimage import io
+
+img = io.imread('data/gaint_panda.bmp')
+imgGray = color.rgb2gray(img)
+M, N = imgGray.shape
+missing_rate = 0.9
+
+sparse_img = imgGray * np.round(np.random.rand(M, N) + 0.5 - missing_rate)
+io.imshow(sparse_img)
+plt.axis('off')
+plt.imsave('gaint_panda_gray_missing_rate_90.png', 
+           sparse_img, cmap = plt.cm.gray)
+plt.show()
+```
+
+```python
+lmbda = 1e-4 * M * N
+gamma = 1 * lmbda
+eta = 100 * lmbda
+tau = 2
+maxiter = 100
+mat_hat = laplacian_conv_2d(imgGray, sparse_img, lmbda, gamma, eta, tau, maxiter)
+
+mat_hat[mat_hat < 0] = 0
+mat_hat[mat_hat > 1] = 1
+io.imshow(mat_hat)
+plt.axis('off')
+plt.imsave('gaint_panda_gray_recovery_90_gamm1_tau{}.png'.format(tau), 
+           mat_hat, cmap = plt.cm.gray)
+plt.show()
+```
+