Robust PCA

When dealing with high-dimensional data sets, we often leverage on the fact that the data has low intrinsic dimensionality in order to alleviate the curse of dimensionality and scale (perhaps it lies in a low-dimensional subspace or lies on a low-dimensional manifold). Principal component analysis

Classical PCA seeks the best rank-k$k$ estimate L$L$ of M$M$ (minimizing ‖M−L‖$\Vert M-L\Vert$ where L$L$ has rank-k$k$). Truncated SVD makes this calculation!

Traditional PCA can handle small noise, but is brittle with respect to grossly corrupted observations– even one grossly corrupt observation can significantly mess up answer.

Robust PCA factors a matrix into the sum of two matrices, M=L+S$M=L+S$, where M$M$ is the original matrix, L$L$ is low-rank, and S$S$ is sparse. This is what we’ll be using for the background removal problem!

• Low-rank means that the matrix has a lot of redundant information– in this case, it’s the background, which is the same in every scene.
  下圖是原矩陣
  M$M$ (4800, 11300)，4800是一幀圖片60x80拉長後的向量長度，11300是幀數。圖中的solid橫線代表background；wave代表行人的moving；column代表 a moment in time。原矩陣 M$M$ 每一列都是某時間點上一幀影象拉長後形成的vector，背景矩陣的rows（圖中的horizontal lines）indicates L$L$ 應當是個low-rank的矩陣；
  
  
• Sparse means that the matrix has mostly zero entries– in this case, see how the picture of the foreground (the people) is mostly empty. (In the case of corrupted data, S$S$ is capturing the corrupted entries. Assume that corrupt data is sparse).

Robust PCA的應用

Face Recognition

Robust PCA Optimization

Source

The general primary component pursuit algorithm from this Robust PCA paper (Candes, Li, Ma, Wright), in the specific form of Matrix Completion via the Inexact ALM Method found in section 3.1 of this paper (Lin, Chen, Ma).

Optimize Algorithms

I. principal component pursuit

Robust PCA can be written:

• ‖⋅‖1$\Vert \cdot {\Vert }_1$ is the L1 norm. Minimizing the L1 norm results in sparse values. For a matrix, the L1 norm is equal to the maximum absolute column norm.

  • Q：這裡為啥用 L1$L_1$ norm呢？

  因為 L1$L_1$ norm會lead更好的稀疏性，而我們需要 S$S$ to be a sparse matrix.

• ‖⋅‖∗$\Vert \cdot {\Vert }_{\ast }$ is the nuclear norm, which is the L1 norm of the singular values. Trying to minimize this results in sparse singular values –> low rank.
  
  

• shrinkage operator Sτ[x]=sgn(x)max(|x|−τ,0)$S_{\tau }[x]=\text{sgn}(x)max(|x|-\tau ,0)$
  取 S$S$ 中的singular value，然後減去 τ$\tau$ ，如果value比 τ$\tau$ 小，就round them to 0. 如果 value 大於 τ$\tau$ 就將兩者的差反號，使得singular value更接近0。總的來說這個shrinkage操作就是使得singular value smaller，從而使它們更加sparse。不要忘了我們的目標是：找到一個sparse的 S$S$。

• singular value thresholder D$D$
  就是先對矩陣做SVD分解，然後通過 S$S$ 的 shrinkage 操作使得singular value變小之後，再reconstruct Dτ(X)=USτ(∑)V∗$D_{\tau }(X)=U{S}_{\tau }(\sum )V^{\ast }$ 得到Low_rank矩陣 L$L$ 的approximation。

  這個操作有點像 truncated SVD，也就是說把singular value比較小（小於 τ$\tau$）的那些vector cut off。

• 第3、4步的意思大概是，因為我們的目標是要把原矩陣分解成一個 Low_rank 矩陣 L$L$ 加上一個 Sparse 矩陣 S$S$ ，所以這裡我們要不斷讓 Low_rank矩陣等於原矩陣 M$M$ 減去 Sparse 矩陣 S$S$，同時讓 S$S$ 儘可能等於 M−L$M-L$。是一個比較典型的交替更新的方法。

• 而 μ−1Yk${\mu }^{-1}{Y}_k$ 是為了 keep track of what you still have to approximate（residual error）.

II. ALM

Section 3.1 of Chen, Lin, Ma contains a faster vartiation of this:

And Section 4 has some very helpful implementation details on how many singular values to calculate (as well as how to choose the parameter values):

Implement principal component pursuit

from scipy import sparse
from sklearn.utils.extmath import randomized_svd
import fbpca

# TOL是收斂的boundary
TOL=1e-9
MAX_ITERS=3

def converged(Z, d_norm):
    err = np.linalg.norm(Z, 'fro') / d_norm
    print('error: ', err)
    return err < TOL

shrinkage operator：把小於 τ$\tau$ 的singular value都置為0，從另一角度來看就是類似在truncating matrix，ignore那些很小的singular value。

def shrink(M, tau):
    S = np.abs(M) - tau
    return np.sign(M) * np.where(S>0, S, 0)

# 寫作pca， 讀作svd，實際上是一個fast svd implementation
def _svd(M, rank): return fbpca.pca(M, k=min(rank, np.min(M.shape)), raw=True)

def norm_op(M): return _svd(M, 1)[1][0]

# D操作
# 這裡的rank是因為fbpca.pca做的是randomized svd，所以需要指定你要保留的rank數
def svd_reconstruct(M, rank, min_sv):
    u, s, v = _svd(M, rank)
    # 這裡面min_sv就是下面提到的sv
    s -= min_sv
    nnz = (s > 0).sum()
    return u[:,:nnz] @ np.diag(s[:nnz]) @ v[:nnz], nnz

principal component pursuit關鍵程式碼實現：

# principal component pursuit是Robust PCA的其中一種algorithm
def pcp(X, maxiter=10, k=10): # refactored
    m, n = X.shape
    trans = m<n
    if trans: X = X.T; m, n = X.shape

    lamda = 1/np.sqrt(m)
    op_norm = norm_op(X)
    Y = np.copy(X) / max(op_norm, np.linalg.norm( X, np.inf) / lamda)
    mu = k*1.25/op_norm; mu_bar = mu * 1e7; rho = k * 1.5

    d_norm = np.linalg.norm(X, 'fro')
    L = np.zeros_like(X); sv = 1

    examples = []

    for i in range(maxiter):
        print("rank sv:", sv)
        X2 = X + Y/mu

        # update estimate of Sparse Matrix by "shrinking/truncating": original - low-rank
        # 對應step 4
        S = shrink(X2 - L, lamda/mu)

        # update estimate of Low-rank Matrix by doing truncated SVD of rank sv & reconstructing.
        # count of singular values > 1/mu is returned as svp
        # 對應step 3
        L, svp = svd_reconstruct(X2 - S, sv, 1/mu)

        '''
        這個sv是個劃重點的東西，代表了你在truncating svd希望保留的rank數目，這個值需要不停地更新
        原始矩陣為400*11300，很明顯我們不可能做full svd，肯定賊慢
        所以我們先甩掉singular value小於1/μ對應的那些vector，
            * 如果我們得到的rank：svp < sv，說明我們想要的value（不小於1/μ的）都已經保留了，
              即我們甩掉的value都是該甩的，則下一輪只將 sv + 1（因為第sv之後的singular value也沒包含太多資訊了）
            * 否則即 svp = sv 說明第sv個singular value之後的value可能還有大於1/μ的，
              我們甩多了，所以下一輪就要多保留一些sv
        看Results中的輸出來check一下
        '''

        # If svp < sv, you are already calculating enough singular values.
        # If not, add 20% (in this case 240) to sv
        sv = svp + (1 if svp < sv else round(0.05*n))

        # residual
        Z = X - L - S
        Y += mu*Z; mu *= rho

        # keep track Time=140 時的visual
        examples.extend([S[140,:], L[140,:]])

        if m > mu_bar: m = mu_bar
        if converged(Z, d_norm): break

    if trans: L=L.T; S=S.T
    return L, S, examples

Results

m, n = M.shape
round(m * .05)

#240

L, S, examples =  pcp(M, maxiter=5, k=10)

rank sv: 1
error:  0.131637937114
rank sv: 241
error:  0.0458515689041
rank sv: 49
error:  0.00588796042886
rank sv: 289
error:  0.000565564416732
rank sv: 529
error:  2.47254909523e-05

sv 初始化為 1，然後加上 n 的 20%，即 240。
然後 svp = 48 < 241，所以 sv = 48 + 1 =49。