Computing the smallest eigenvalue of a symmetric matrix

In this example we compute the smallest eigenvalue of a symmetric matrix $A \in ℝ^{n×n}$. Let $λ_1 ≥ ... ≥ λ_n$ be the eigenvalues of $A$. Then the smallest eigenvalue $λ_n$ can be computed by solving the following optimization problem:

$$ \operatorname*{arg\,min}_{x∈ℝ^n, x ≠ 0} \frac{x^{\mathrm{T}}Ax}{x^{\mathrm{T}}x}.$$

This problem can be reformulated as

$$\operatorname*{arg\,min}_{x∈ℝ^n, \lVert x\rVert_2^2 = 1} x^{\mathrm{T}}Ax$$

and thus becomes an optimization problem on the sphere $𝕊^{n-1} = \{ x∈ℝ^n : \lVert x\rVert _2^2 = x^{\mathrm{T}}x = 1 \}.$

The cost function and its gradient in $ℝ^n$ can be written as: $f(x) = x^{\mathrm{T}}Ax$ and $\nabla f(x) = 2Ax$.

For the optimization on the manifold $𝕊^{n-1}$ we also need the Riemannian gradient of $f$. It can be calculated by applying the orthogonal projection onto the tangent space of $x$ $\operatorname{proj}_x u = (I-xx^{\mathrm{T}})u$ to $∇f(x)$:

$$ \operatorname{grad}f(x) = \operatorname{proj}_x(∇f(x)) = 2 (I-xx^{\mathrm{T}}) Ax.$$

First, for illustration purposes, we create a random matrix $A$ whose smallest eigenvalue is then computed.

using Random, LinearAlgebra, PlutoUI
Random.seed!(42);
n = 10;
a = randn(n, n);
A = a * a';

For the solution of the optimization problem we want to use the Manopt toolbox. For this purpose we now create our manifold $𝕊^{n-1}$ and a starting vector $x_0 ∈ 𝕊^{n-1}$.

We also implement the cost function f and its Riemannian gradient gradf using the Sphere from Manifolds.jl.

using Manifolds
x0 = ones(n) ./ sqrt(n)
10-element Vector{Float64}:
 0.31622776601683794
 0.31622776601683794
 0.31622776601683794
 0.31622776601683794
 0.31622776601683794
 0.31622776601683794
 0.31622776601683794
 0.31622776601683794
 0.31622776601683794
 0.31622776601683794

Here we can first check that its norm (0.9999999999999999) is numerically one or in other words is_point(M,x0) yields true

f(M, x) = x' * A * x
f (generic function with 1 method)
gradf(M, x) = 2.0 * (I(n) - x * x') * A * x
gradf (generic function with 1 method)
M = Sphere(n - 1)
Sphere(9, ℝ)

For optimization we use a gradient descent method with Armijo linesearch. The selected linesearch ensures a sufficient descent.

using Manopt
with_terminal() do
    global xopt = gradient_descent(
        M,
        f,
        gradf,
        x0;
        stepsize=ArmijoLinesearch(M; contraction_factor=0.99, sufficient_decrease=0.6),
        debug=[:Iteration, :Cost, ", ", DebugGradientNorm(), 50, :Stop, "\n"],
    )
end
Initial F(x): 12.524120, 
# 50    F(x): 0.067004, |gradF(x)|:0.020201948409546262
# 100   F(x): 0.066849, |gradF(x)|:9.097563382427176e-5
The algorithm reached approximately critical point after 132 iterations; the gradient norm (9.492360108863904e-9) is less than 1.0e-8.

Now we check whether the result of the optimization is reasonable. For this purpose, we use the objective function value of the optimizer and, for comparison, a function which provides the minimum eigenvalue of a matrix.

Furthermore, the Euclidean norm of the optimizer is calculated to check whether the point lies on the sphere.

We first compute as the “ground truth”

t = minimum(eigvals(A))
0.06684902616321216

Compared to the cost at xopt, i.e.

v = f(M, xopt)
0.0668490261631732

which is a good approximation of the original eigenvalue:

v - t
-3.895495037653518e-14

And we can also check whether xopt is a point on the sphere:

is_point(M, xopt)
true