1. 程式人生 > >Robust PCA——Inexect ALM(含python實現)

Robust PCA——Inexect ALM(含python實現)

原文:https://blog.csdn.net/u010510350/article/details/77885553
前兩篇部落格已經介紹了Robust PCARPCA的優化,接下來用Robust PCA實現背景建模。背景建模就是將攝像機獲取的場景分離出前景和背景,以獲取場景中的動態目標。傳統方法基本思路:首先通過學習一段訓練影象序列提取出該視訊的背景特徵,來建立數學模型以便描述其背景,然後用該背景模型對需要檢測的視訊序列進行處理( 如採用背景相減法) ,最後提取出當前影象中與背景模型中性質不同的畫素點,即為影象的動態目標。但由於視訊場景中可能存在噪聲汙染,所以傳統的背景建模方式就沒法滿足建模要求。Robust PCA,又稱稀疏與低秩矩陣分解,能夠從受強噪聲汙染或部分缺失的高維度觀測樣本中發現其低維特徵空間,有效地恢復觀測樣本的低維子空間,並恢復受損的觀測樣本。

1.背景建模Robust PCA

,mVi,n,nD=[v1,v2,...,vn]Rmn(i=1,2,...,n)RobustPCA:D,(A)(E):對於某類觀測的視訊影象,將每一幀影象表示為m維向量Vi,若該視訊共包含n幀影象序列,那麼該觀測視訊就可以用n個向量組成的資料矩陣D=[v1,v2,...,vn]∈Rm∗n(i=1,2,...,n)來表示。視訊背景建模的RobustPCA可描述為:對資料矩陣D進行分解,恢復出具有極大相似性的背景部分(低秩矩陣A)和分佈範圍很小的運動目標或前景部分(稀疏矩陣E):

這裡寫圖片描述
這樣就形成了yongqianggao中的凸優化問題,這樣就可以通過上篇部落格介紹的優化方法進行求解了,計算出A和E就可以分離出背景和前景(運動目標)了。

2.Inexect ALM求解A和E

上篇部落格已經介紹了Robust PCA的Inexect ALM優化演算法,在這裡就不進行贅述,直接上馬毅2009年”The Augmented Lagrange Multiplier Method for Exact Recovery of Corrupted Low-Rank Matrices“論文中演算法截圖。
這裡寫圖片描述

3.Inexect ALM演算法實現

下面是IALM演算法的Python實現,具體可參考我的分享的

程式碼資源

import numpy as np 
from numpy.linalg import norm, svd

def inexact_augmented_lagrange_multiplier(X, lmbda = 0.01, tol = 1e-7, maxIter = 1000):
    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[1]
    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(0.05 * n), n])

        Aupdate = np.dot(np.dot(U[:, :svp], np.diag(S[:svp] - 1 / mu)), V[:svp, :])
        A = Aupdate
        E = Eupdate
        print itr
        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
    print("IALM Finished at iteration %d" % (itr))
    return A, E