欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

Robust PCA

程序员文章站 2022-05-05 08:49:27
...

拉格朗日乘子解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)
相关标签: 预处理