Robust PCA

拉格朗日乘子解Robust PCA以及Python實作

————————————————

版權聲明：本文為CSDN部落客「masonwang_513」的原創文章，遵循CC 4.0 BY-SA版權協定，轉載請附上原文出處連結及本聲明。

原文連結：https://blog.csdn.net/reform513/article/details/79539511

Implementation of Python

import numpy as np
from numpy.linalg import norm, svd
 
def inexact_augmented_lagrange_multiplier(X, lmbda=.01, tol=1e-3,
                                          maxiter=100, verbose=True):
    """
    Inexact Augmented Lagrange Multiplier
    """
    Y = X
    norm_two = norm(Y.ravel(), 2)
    norm_inf = norm(Y.ravel(), np.inf) / lmbda
    dual_norm = np.max([norm_two, norm_inf])
    Y = Y / dual_norm
    A = np.zeros(Y.shape)
    E = np.zeros(Y.shape)
    dnorm = norm(X, 'fro')
    mu = 1.25 / norm_two
    rho = 1.5
    sv = 10.
    n = Y.shape[0]
    itr = 0
    while True:
        Eraw = X - A + (1 / mu) * Y
        Eupdate = np.maximum(Eraw - lmbda / mu, 0) + np.minimum(Eraw + lmbda / mu, 0)
        U, S, V = svd(X - Eupdate + (1 / mu) * Y, full_matrices=False)
        svp = (S > 1 / mu).shape[0]
        if svp < sv:
            sv = np.min([svp + 1, n])
        else:
            sv = np.min([svp + round(.05 * n), n])
        Aupdate = np.dot(np.dot(U[:, :svp], np.diag(S[:svp] - 1 / mu)), V[:svp, :])
        A = Aupdate
        E = Eupdate
        Z = X - A - E
        Y = Y + mu * Z
        mu = np.min([mu * rho, mu * 1e7])
        itr += 1
        if ((norm(Z, 'fro') / dnorm) < tol) or (itr >= maxiter):
            break
    if verbose:
        print("Finished at iteration %d" % (itr))  
    return A, E
           

num, d1, d2 = 500, 180, 320
dt = np.load('test\\image_data.npz')['data']
data = np.empty((num, d1, d2), dtype=np.uint8)
for i in range(num):
    data[i] = cv2.resize(dt[i], dsize=(d2, d1))
X = data.reshape(num, -1).T / 255.0
print(X.shape)
A, E = inexact_augmented_lagrange_multiplier(X)
 
A = (A.T.reshape(num, d1, d2) * 255).astype('uint8')
E = (E.T.reshape(num, d1, d2) * 255).astype('uint8')
for i in range(data.shape[0]):
    cv2.imshow('Original', cv2.resize(data[i], dsize=(d2, d1)))
    cv2.imshow('LowRank',  cv2.resize(A[i], dsize=(d2, d1)))
    cv2.imshow('Sparse',   cv2.resize(E[i], dsize=(d2, d1)))
    cv2.waitKey(0)