I recently spent some time working on a simple linear algebra problem - decompose a matrix $ M $ into a low-rank component $ L $ and a sparse component $ S $. The algorithm I used was very trivial to implement (and parallelize using map-reduce).
In this post, I will implement this very simple algorithm, explain the objective function and demonstrate its (amazing) effectiveness on a surveillance-camera dataset.
Given a matrix $ M $, we want to decompose it into $ L + S $ where $ L $ is a dense, low-rank matrix and $ S $ is a sparse matrix. The typical approach to solve these problems involves solving an optimization problem involving the nuclear norm.
Surprisingly, the PCP estimate fully recovers both these matrices. Its objective is $ || L_{0} ||_{*} + \lambda||S||_{1} $ subject to $ L + S = M $. Essentially, the $ l1 $ norm enforces sparsity in S and the constraint is intuitive. This works because:
If the singular vectors are reasonably spread, they can be recovered with very high accuracy. This is clearly seen since a whole host of algorithms we use on a daily basis make use of this. The recovery is also not perturbed too much by slight corruptions (which is related to the sparse component). Thus, recovering the singular values while enforcing sparsity in the corruption matrix (the sparse matrix) is a sufficient condition for recovering the original matrix.
Also, the lambda component itself needs no tuning (!!!!). So no cross-validation component exists. $ \lambda $ is just set to $ \frac{1}{\sqrt{n}} $ where $ n $ is the dimension of the matrix.
The augmented Lagrangian problem works with the following:
$ l(L, S, Y) = ||L||_{*} + \lambda||S||_{1} + \langle Y, M - L - S \rangle + \frac{\mu}{2} ||M - L - S||_{F}^{2} $
where:
- The $ ||L||_{*} $ is the nuclear norm and this ensures that the matrix $ L $ is low-rank
- The $ ||S||_{1} $ is sparsity inducing $ l1 $ norm.
- The other two expressions drive the value $ L + S $ towards that of $ M $.
The iteration steps are:
init: S0 = Y0 = 0, mu > 0
while not converged {
Li+1 = SVD_thresh(M - Si - (1/mu)Yi, mu);
Si+1 = thresh(M - Li+1 + (1/mu)Yi, lambda * mu);
Yi+1 = Yi + mu(M - Li+1 - Si+1);
}
return L, Swhere SVD_thresh reconstructs a matrix by zero’ing out singular values that fall below a certain threshold and thresh zeros out matrix elements that fall below the specified threshold.
This simple algorithm accurately recovers the low-rank and sparse matrix components!
Background Detection
A surveillance video camera emits a sequence of frames. Now, consider a matrix where one single row (or an individual column) is created using all the pixels of a frame (so a 160 x 130 video becomes a 20800 dimension vector). Since the background component is common across frames, the low-rank component is just a rank–1 matrix (the only basis vector is the background itself). The sparse component consists of the moving parts (people, escalators, walkways and so on).
Pasted below is the python version of the algorithm above (this is a port from the Matlab version shared by the authors here and uses the excellent PROPACK library to compute the partial SVD):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 |
#!/usr/bin/env python import sys, os import argparse import numpy as np from pypropack import svdp; from scipy.io import loadmat MAX_ITERS = 1000 TOL = 1.0e-7 def process_args(): parser = argparse.ArgumentParser() parser.add_argument( '--mat-file', metavar = 'mat_file', help = 'Location to .mat file' ) parser.add_argument( '--test', dest = 'test', default = False, action = 'store_true' ) parser.add_argument( ) return parser.parse_args() def choose_svd(n, sv): return True def converged(Z, d_norm): err = np.linalg.norm(Z, 'fro') / d_norm print 'ERR', err return err < TOL def pcp(X): m, n = X.shape # Set params lamda = 1./ np.sqrt(m); # Initialize Y = X; u, s, v = svdp(Y, k=1, which='L'); norm_two = s[0] norm_inf = np.linalg.norm( Y[:], np.inf) / lamda; dual_norm = max(norm_two, norm_inf); Y = Y / dual_norm; A_hat = np.zeros((m, n)); E_hat = np.zeros((m, n)); mu = 1.25/norm_two mu_bar = mu * 1e7 rho = 1.5 d_norm = np.linalg.norm(X, 'fro'); num_iters = 0; total_svd = 0; stopCriterion = 1; sv = 10; while True: num_iters += 1; temp_T = X - A_hat + (1/mu)*Y; E_hat = np.maximum(temp_T - lamda/mu, 0); E_hat = E_hat + np.minimum(temp_T + lamda/mu, 0); u, s, v = svdp(X - E_hat + (1/mu)*Y, sv, which = 'L'); diagS = np.diag(s); svp = len(np.where(s > 1/mu)) if svp < sv: sv = min(svp + 1, n); else: sv = min(svp + round(0.05*n), n); A_hat = np.dot( np.dot( u[:,0:svp], np.diag(s[0:svp] - 1/mu) ), v[0:svp,:] ) total_svd = total_svd + 1; Z = X - A_hat - E_hat; Y = Y + mu*Z; mu = min(mu*rho, mu_bar); if converged(Z, d_norm) or num_iters >= MAX_ITERS: return A_hat, E_hat if __name__ == '__main__': args = process_args() # Load Data if not args.test: data = loadmat(args.mat_file) A_hat, E_hat = pcp(data['X']) print A_hat, E_hat else: data = (10*np.ones((10, 10))) + (-5 * np.eye(10)) print data A_hat, E_hat = pcp(data) np.save('low_rank', A_hat) np.save('sparse', E_hat) |
Next, let us deploy this on a dataset. I obtained a video from a surveillance camera (a few datasets are available at this webpage. I will use the Escalator dataset from that page).
The following script builds a numpy matrix from the bitmap image sequence:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
''' Convert a bitmap sequence to matrix. ''' import os import numpy as np import argparse from PIL import Image def bitmap_to_mat(bitmap_seq): matrix = [] for bitmap_file in bitmap_seq: im = Image.open(bitmap_file).convert('L') # convert to grayscale pixels = list(im.getdata()) matrix.append(pixels) return np.array(matrix) if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('--bmp-dir', dest = 'bmp_dir') parser.add_argument('--bmp-out', dest = 'bmp_out') parsed = parser.parse_args() bmp_seq = map( lambda s: os.path.join(parsed.bmp_dir, s), os.listdir(parsed.bmp_dir) ) res = bitmap_to_mat(bmp_seq) np.save(parsed.bmp_out, res) |
Next, let us run the decomposition and see what the convergence looks like:
➜ robust_pcp git:(master) python robust_pcp.py --input-file bmp_matrix.npy
ERR 0.799466063211
WARNING: Maximum dimension of Krylov subspace exceeded prior to convergence. Try increasing KMAX.
neig = 1
ERR 0.533151636032
WARNING: Maximum dimension of Krylov subspace exceeded prior to convergence. Try increasing KMAX.
neig = 1
ERR 0.580203519826
ERR 0.399795997647
WARNING: Maximum dimension of Krylov subspace exceeded prior to convergence. Try increasing KMAX.
neig = 1
ERR 0.526129524054
WARNING: Maximum dimension of Krylov subspace exceeded prior to convergence. Try increasing KMAX.
neig = 1
ERR 0.576362875779
ERR 0.438452861276
ERR 0.137666099277
WARNING: Maximum dimension of Krylov subspace exceeded prior to convergence. Try increasing KMAX.
neig = 1
ERR 0.526604907584
WARNING: Maximum dimension of Krylov subspace exceeded prior to convergence. Try increasing KMAX.
neig = 1
ERR 0.57572793864
ERR 0.45163686726
ERR 0.111636379698
ERR 0.0718803322749
ERR 0.048067444228
WARNING: Maximum dimension of Krylov subspace exceeded prior to convergence. Try increasing KMAX.
neig = 1
ERR 0.536409101582
ERR 0.373332905257
ERR 0.0414766101901
ERR 0.0264232708867
WARNING: Maximum dimension of Krylov subspace exceeded prior to convergence. Try increasing KMAX.
neig = 1
ERR 0.509857058194
ERR 0.341770886808
ERR 0.00828382155722
ERR 0.00369347648645
ERR 0.00224842770609
ERR 0.00133359512154
ERR 0.00076422430761
ERR 0.00043106684577
ERR 0.000242042207018
ERR 0.000134750434297
ERR 7.47073010806e-05
ERR 4.12146122072e-05
ERR 2.26314422366e-05
ERR 1.23645958329e-05
ERR 6.76884094823e-06
ERR 3.71124761483e-06
ERR 2.01918022651e-06
ERR 1.10232817344e-06
ERR 6.04773368406e-07
ERR 3.27599538767e-07
ERR 1.79967651854e-07
ERR 9.8752813134e-08So in 40 iterations, the low-rank and sparse matrices are completely recovered.
Now, let us take a look at what the decomposition looks like. The following piece of code take one of those matrices we built (from the bitmaps) and animates them (i.e. iterates through the sequence):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 |
''' Animate supplied npy files ''' import argparse import numpy as np import time import Tkinter from PIL import Image, ImageTk def animate(matrices, w, h): mats = [np.load(x) for x in matrices] tk_win = Tkinter.Toplevel() tk_win.title('Escalator') canvas = Tkinter.Canvas(tk_win, width=7*w, height=7*h) canvas.pack() tk_ims = [None for _ in mats] for i, row in enumerate(mats[0]): ims = [Image.new('L', (w, h)) for _ in mats] for j, im in enumerate(ims): im.putdata(map(float, list(mats[j][i]))) tk_ims[j] = ImageTk.PhotoImage(im) canvas.create_image((j * w) + 200, h, image = tk_ims[j]) canvas.update() if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument( 'width', type = int, help = 'frame width' ) parser.add_argument( 'height', type = int, help = 'frame width' ) parser.add_argument( '--npy-files', dest = 'npy_files', nargs = '+', help = 'numpy matrices' ) parsed = parser.parse_args() animate(parsed.npy_files, parsed.width, parsed.height) |
And here’s a frame from the resulting animation. The first frame is the original image sequence (in grayscale), the second is the low-rank component (the background) and the sparse component consists of the people:
The algorithm is thus super-powerful and easy to implement.
The full source code is available in this github repository.
References: