# Robust PCA (better Mathematica style)

I’ve parsed (almost verbatim) Python RPCA implementation to WM. Can it be rewritten using a better WM style? In particular, I’m not happy with `Clip[..., {0, Infinity}]` and `While[]` loop (can be replaced with `Do[]`).

``ClearAll[shrink] ; shrink[matrix_, tau_] := Sign[matrix]*Clip[Abs[matrix] - tau, {0, Infinity}]  ClearAll[threshold] ; threshold[matrix_, tau_] := Block[     {u, s, v},     {u, s, v} = SingularValueDecomposition[matrix] ;     Dot[u, Dot[shrink[s, tau], Transpose[v]]] ] ;  ClearAll[rpca] ; rpca[matrix_?MatrixQ, mu_Real, lambda_Real, tolerance_Real, limit_Integer] := Block[     {inverse, count, error, sk, yk, lk},     inverse = 1.0/mu ;      count = 0 ;     sk = yk = lk = ConstantArray[0.0, Dimensions[matrix]] ;     While[         count < limit,         lk = threshold[matrix - sk + inverse*yk, inverse] ;         sk = shrink[matrix - lk + inverse*yk, inverse*lambda] ;         error = matrix - lk - sk ;         yk = yk + mu*error ;         error = Norm[error, "Frobenius"] ;         count++ ;         If[error < tolerance, Break[]] ;     ] ;     {lk, sk, {count, error}} ] ; ``

Example:

``(* https://github.com/dganguli/robust-pca *) (* "12.1.1 for Linux x86 (64-bit) (June 19, 2020)" *)  (* generate test matrix *)  n = 100 ; num$$groups = 3 ; num$$values = 40 ; matrix = N[ConstantArray[Flatten[Transpose[ConstantArray[10*Range[num$$groups], num$$values]]], n]] ; {n, m} = Dimensions[matrix]  (* set selected elements to zero *)  SeedRandom[1] ; ln = RandomInteger[{1, n}, 20] ; lm = RandomInteger[{1, m}, 20] ; Table[matrix[[ln[[i]], lm[[i]]]] = 0.0, {i, 1, 20}] ; matrix = Developer`ToPackedArray[matrix] ;  (* -- python zeros ln = [81, 15, 1, 68, 4, 66, 24, 98, 69, 75, 16, 25, 5, 91, 84, 71, 2, 31, 49, 26] lm = [45, 74, 107, 70, 57, 48, 29, 69, 27, 69, 11, 87, 77, 44, 34, 45, 87, 94, 19, 39] for x, y in zip(ln, lm):     D[x, y] = 0 *)  (* set parameters *) mu = 1/4*1/Norm[matrix, 1]*Apply[Times, Dimensions[matrix]] ; lambda = 1/Sqrt[N[Max[Dimensions[matrix]]]] ; tolerance = 10.0^-7*Norm[matrix, "Frobenius"] ; limit = 1000 ;  (* rpca *) result = rpca[matrix, mu, lambda, tolerance, limit] ;  (* # of iterations and error *) Last[result]  (* low rank *) +Table[result[[1]][[ln[[i]], lm[[i]]]], {i, 1, 20}]  (* sparse *) -Table[result[[2]][[ln[[i]], lm[[i]]]], {i, 1, 20}]  (* {100, 120} *) (* {39, 0.000167548} *) (* {20., 20., 30., 20., 20., 20., 10., 20., 10., 20., 10., 23.0449, 20., 20., 10., 20., 23.0449, 30., 10., 10.} *) (* {20., 20., 30., 20., 20., 20., 10., 20., 10., 20., 10., 23.0449, 20., 20., 10., 20., 23.0449, 30., 10., 10.} *) ``