# Computing a Robust PCA

For a given matrix $X ∈ ℝ^{p×n}$ whose columns represent points in $ℝ^p$, a matrix $U ∈ ℝ^{p×d}$ is computed for a given dimension $d < n$: $U$ represents an ONB of $ℝ^{p× d}$ such that the column space of $U$ approximates the points (columns of $X$) $X_i$ as well as possible.

We compute $U$ as a minimizer over the Grassmann manifold of the cost function:

$$\begin{split} f(U) & = \frac{1}{n}\sum_{i=1}^{n}{\operatorname{dist}(X_i, \operatorname{span}(U))}\\ & = \frac{1}{n} \sum_{i=1}^{n}\lVert UU^TX_i - X_i\rVert \end{split}$$

The output cost represents the average distance achieved with the returned $U$, an orthonormal basis (or a point on the Stiefel manifold) representing the subspace (a point on the Grassmann manifold). Notice that norms are not squared, so we have a robust cost function. This means that $f$ is nonsmooth, therefore we regularize with a pseudo-Huber loss function of smoothing parameter $ϵ$.

$$f_ϵ(U) = \frac{1}{n} \sum_{i=1}^n{ℓ_ϵ(\lVert UU^{\mathrm{T}}X_i - X_i\rVert)},$$

Where $ℓ_ϵ(x) = \sqrt{x^2 + ϵ^2} - ϵ$.

The smoothing parameter is iteratively reduced (with warm starts).

First, we generate random data for illustration purposes:

using Random
begin
samples = 40
outliers = 15
Random.seed!(42)
X = randn(2, 1) * [1:samples;]' + 0.05 * randn(2, samples) .* [1:samples 1:samples]'
# Outliers:
P = shuffle(1:size(X, 2))'
X[:, P[1:outliers]] = 30 * randn(2, outliers)
d = 1
end;

We use the Manopt toolbox to optimize the regularized cost function over the Grassmann manifold.

To do this, we first need to define the problem structure.

using Manifolds
p, n = size(X);
M = Grassmann(p, d);

For the initial matrix $U_0$ we use classical PCA via singular value decomposition. Thus, we use the first $d$ left singular vectors.

Then, we compute an optimum of the cost function over the Grassmann manifold. We use a trust-region method which is implemented in Manopt.jl.

For this, the objective function robustpca_cost and the Riemannian gradient grad_pca_cost must first be implemented as functions. Since we Huber regularize them, the functions also have the ϵ as a parameter. To obtain the Riemannian gradient we first compute the Euclidian gradient. Afterwards it is projected onto the tangent space by using the orthogonal projection $I-U*U^T$.

The trust-region method also requires the Hessian Matrix. By using ApproxHessianFiniteDifference we get an approximation of the Hessian Matrix.

We run the procedure several times, where the smoothing parameter $ϵ$ is reduced iteratively.

using LinearAlgebra, Manopt, Statistics
begin
ϵ = 1.0
iterations = 6
reduction = 0.5
end;
U, S, V = svd(X);
U0 = U[:, 1:d];
U0
2×1 Matrix{Float64}:
-0.7494248652139397
0.6620893983436593
robustpca_cost(M, U0, 0.0)
10.218786443576352
function grad_pca_cost(M, U, ϵ)
UtX = U' * X
vecs = U * UtX - X
sqnrms = sum(vecs .^ 2; dims=1)

G = zeros(p, d)
for i in 1:n
G = G + (1 / (sqrt(sqnrms[i] + ϵ^2))) * vecs[:, i] * UtX[:, i]'
end
G = 1 / n * G
return (I - U * U') * G
end
grad_pca_cost (generic function with 1 method)
function robustpca_cost(M, U, ϵ)
vecs = U * U' * X - X
sqnrms = sum(vecs .^ 2; dims=1)
dim = size(sqnrms)
vals = sqrt.(sqnrms + ϵ^2 * ones(dim)) - ϵ * ones(dim)
return mean(vals)
end
robustpca_cost (generic function with 1 method)
begin
Ui = copy(M, U0)
ϵi = ϵ
for i in 1:iterations
global Ui = trust_regions(
M,
(M, p) -> robustpca_cost(M, p, ϵi),
(M, p) -> grad_pca_cost(M, p, ϵi),
ApproxHessianFiniteDifference(
M,
Ui,
(M, p) -> grad_pca_cost(M, p, ϵi);
vector_transport_method=ProjectionTransport(),
retraction_method=PolarRetraction(),
),
Ui;
(project!)=project!,
)
global ϵi = ϵi * reduction
end
end;
Ui
2×1 Matrix{Float64}:
-0.6795640993446082
0.7336161359198359
robustpca_cost(M, Ui, 0.0)
9.412961993577344

Finally, the results are presented graphically. The data points are visualized in a scatter plot. The result of the robust PCA and (for comparison) the standard SVD solution are plotted as straight lines.

using Plots
fig = plot(X[1, :], X[2, :]; seriestype=:scatter, label="Data points");
plot!(
fig,
Ui[1] * [-1, 1] * 100,
Ui[2] * [-1, 1] * 100;
linecolor=:red,
linewidth=2,
label="Robust PCA",
);
plot!(
fig,
U0[1] * [-1, 1] * 100,
U0[2] * [-1, 1] * 100;
xlims=1.1 * [minimum(X[1, :]), maximum(X[1, :])],
ylims=1.1 * [minimum(X[2, :]), maximum(X[2, :])],
linewidth=2,
linecolor=:black,
label="Standard SVD",
)