Sparse principal components

Sparse principal component analysis by Fantope projection and selection.

maximize  ⟨S, X⟩ - λ ‖X‖₁   subject to  X ∈ Fᵏ = {X : 0 ⪯ X ⪯ I, tr(X) = k}

The Fantope Fᵏ is the convex hull of the projections onto k-dimensional subspaces, so without the penalty the maximizer is the projection onto the leading k eigenvectors of the sample covariance S. The ℓ1 penalty makes the estimate of the subspace sparse (Vu, Cho, Lei and Rohe, 2013).

The indicator of Fᵏ is Spectral(Constrained(IndBox(0, 1), IndHyperplane(1, k))): its projection projects the eigenvalues onto the capped simplex. ADMM splits the problem as h(X) + g(Y) with X = Y, for h = Tilt(indicator of Fᵏ, -S), whose proximal operator is the projection of X + S/ρ, and g = λ‖⋅‖₁.

using ProximalBase, LinearAlgebra, Random

rng = Xoshiro(17)

Random.Xoshiro(0x86e4c19d9405fc84, 0x43b6af16d34e2a65, 0x5914b730eb37944e, 0x02667705420b1d4a, 0x203da3fb74542e3f)

Data: a spiked covariance with a sparse leading eigenvector

p, n, k = 50, 200, 1
v = zeros(p)
v[1:5] .= [0.6, 0.5, 0.4, 0.35, 0.3]
v ./= norm(v)
Σ = I + 2.0 * v * v'                             # the leading eigenvector is v
X = randn(rng, n, p) * cholesky(Σ).U              # rows from N(0, Σ)
S = Symmetric(X' * X / n)

50×50 LinearAlgebra.Symmetric{Float64, Matrix{Float64}}:
  1.67348     0.717647    0.515158    …  -0.0318799   -0.0200226
  0.717647    1.58155     0.597751       -0.103707    -0.02732
  0.515158    0.597751    1.23933        -0.0576458   -0.0601376
  0.513567    0.35081     0.296302        0.10167     -0.0235139
  0.341465    0.329261    0.168689        0.0508923    0.0222485
  0.123804   -0.0507365   0.00977792  …  -0.109161     0.109117
 -0.102047   -0.0807422  -0.0474026      -0.00839483   0.0561052
 -0.0485354   0.122492    0.0557524      -0.00560918  -0.104943
  0.0591334  -0.0670024  -0.0298017       0.0506705    0.0493601
 -0.0955585  -0.11947    -0.0292286       0.0212479   -0.0887009
  ⋮                                   ⋱               
 -0.217358   -0.0265814  -0.107644        0.137814    -0.0116244
  0.0984075   0.106963   -0.0833812      -0.0812493   -0.0449529
 -0.0170181  -0.0804936  -0.0909464      -0.0518178    0.035773
  0.013062    0.0420517  -0.01591         0.0293147   -0.169642
  0.010331    0.0217644   0.0645598   …  -0.0540304   -0.147969
 -0.0985846  -0.0307143   0.0188641       0.0217222    0.0148886
 -0.0860568  -0.0527603   0.0492549       0.120862    -0.050211
 -0.0318799  -0.103707   -0.0576458       1.02265      0.0466268
 -0.0200226  -0.02732    -0.0601376       0.0466268    0.901578

Plain PCA

u = eigen(S).vectors[:, end]
println("PCA: ", count(x -> abs(x) > 1e-3, u), " nonzero loadings, |⟨u, v⟩| = ", round(abs(dot(u, v)), digits = 3))

PCA: 50 nonzero loadings, |⟨u, v⟩| = 0.934

Fantope projection and selection by ADMM

λ = 0.15
fantope = Spectral(Constrained(IndBox(0, 1), IndHyperplane(1, k)))
h = Tilt(fantope, -Matrix(S))                    # -⟨S, X⟩ on the Fantope
g = NormL1(λ)
ρ = 1.0
Xk, Y, U = zeros(p, p), zeros(p, p), zeros(p, p)
iters = 0
for it in 1:2000
    global iters = it
    prox!(Xk, h, Y .- U, 1 / ρ)                  # the projection of Y - U + S/ρ onto Fᵏ
    Yold = copy(Y)
    prox!(Y, g, Xk .+ U, 1 / ρ)                  # soft thresholding
    U .+= Xk .- Y
    r, s = norm(Xk - Y), ρ * norm(Y - Yold)      # the primal and dual residuals
    max(r, s) < 1e-9 && break
end
println("ADMM: ", iters, " iterations, ‖X - Y‖ = ", norm(Xk - Y))

μ = eigvals(Symmetric(Xk))
@assert all(-1e-9 .<= μ .<= 1 + 1e-9) && tr(Xk) ≈ k        # X is in the Fantope
@assert norm(Xk - Y) < 1e-8

ADMM: 103 iterations, ‖X - Y‖ = 9.490710993897126e-10

the leading eigenvector of the estimate: sparse, and nearer to v than PCA's

w = eigen(Symmetric(Y)).vectors[:, end]
support = findall(x -> abs(x) > 1e-6, w)
println("Fantope projection and selection: loadings on ", support, ", |⟨w, v⟩| = ", round(abs(dot(w, v)), digits = 3))
@assert support == 1:5
@assert abs(dot(w, v)) > abs(dot(u, v))

Fantope projection and selection: loadings on [1, 2, 3, 4, 5], |⟨w, v⟩| = 0.985

the optimality conditions, through the scaled dual variable U: ρU is a subgradient of λ‖⋅‖₁ at Y, and S - ρU is normal to Fᵏ at X, so that the projection of X + S - ρU onto Fᵏ is X

@assert all(abs(ρ * U[i] - λ * sign(Y[i])) < 1e-6 for i in eachindex(Y) if !iszero(Y[i]))
@assert all(abs(ρ * U[i]) <= λ + 1e-6 for i in eachindex(Y) if iszero(Y[i]))
@assert norm(prox(fantope, Xk .+ (Matrix(S) .- ρ .* U)) - Xk) < 1e-6