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): 5.602170, 
# 50    F(x): 0.024540, |gradF(x)|:0.008731640502240161
# 100   F(x): 0.024536, |gradF(x)|:4.958520050014457e-6
# 150   F(x): 0.024536, |gradF(x)|:2.0932355310337282e-7
# 200   F(x): 0.024536, |gradF(x)|:2.0932187085208932e-7
The algorithm reached its maximal number of iterations (200).

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.024536188677448063

Compared to the cost at xopt, i.e.

v = f(M, xopt)

0.0245361886774568

which is a good approximation of the original eigenvalue:

v - t

8.7360674250192e-15

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

is_point(M, xopt)

true