Is there a manifold for Hermitian positive definite matrices?

[Anne Greenbaum]

Given an nxn matrix M, I want to solve the following problem:

min_{X (complex) Hermitian positive definite} || X^{-1} M X || s.t. || X^{-1} || ||X|| <= 2.

Here || || is the operator 2-norm.  Haven't used manopt before.  Any hints about how

to do it appreciated.

[Nicolas Boumal]

Hello,

This is not a typical use case for manopt because the constraint ||X^{-1}|| ||X|| <= 2 does not define a smooth manifold.

Nevertheless, when this happens, one can often still "fit" the problem into a smooth manifold, e.g., by smooth parameterization.

Thinking along those lines, here is one idea that comes to mind: we could parameterize X as X = UDU* with U a unitary matrix (available in manopt as unitaryfactory(n)) and D a diagonal matrix with special diagonal entries.

Here, we would want the diagonal entries of D to be real numbers in the interval [1, 2]. That's because the cost function and constraint are invariant under scaling (so we can fix the lower-bound arbitrarily; here: 1), and the constraints require the ratio between min and max eigenvalue to be at most 2.

The constraint "D_ii is in [1, 2]" is also not a smooth manifold constraint, but it is easily parameterized as D_ii = 1.5 + .5*sin(e), where e is an unconstrained real number. We need n such numbers (for the diagonal of D), so e should really be a vector in R^n. That's available in manopt as euclideanfactory(n, 1).

We can make a product manifold for this:

n = 5; % choose a size

elems.U = unitaryfactory(n);

elems.e = euclideanfactory(n, 1);

manifold = productmanifold(elems);

Then, a point X on the product manifold is a structure with fields X.U (unitary matrix) and X.e (free vector) such that the matrix is:

X_matrix = X.U * diag(1.5 + .5*sin(X.e)) * X.U' ;

For example, this can be checked with:

X = manifold.rand(); % generate a random point on the product manifold

X_matrix = X.U * diag(1.5 + .5*sin(X.e)) * X.U' ;

eig(X_matrix) % by design, all in [1, 2]

As desired, X_matrix is clearly Hermitian positive definite, and its eigenvalues are between 1 and 2, for a condition number of at most 2.

Now, for the cost function f(X) = || X^{-1} M X ||, there is another difficulty because it, too, is not smooth, yet the best algorithms in manopt are designed for smooth optimization. Still, we can try a few things.

Here, the cost function in terms of U and e (from our new parameterization) will be g(U, e) = f(U * diag(1.5 + .5*sin(e)) * U') = || diag((1.5 + .5*sin(e)).^(-1)) * U'*M*U * diag(1.5 + .5*sin(e)) ||  (using the fact that the operator norm is unitarily invariant).

From here, I would start by implementing that cost function to gain some insight, and see from here how much more we want to try that approach.

Here is a complete script -- as a fair warning: this is unlikely to be satisfactory, because the nonsmooth cost function does not allow access to the best solvers.

clear; clc;

n = 5;

M = randn(n) + 1i*randn(n); % general n-by-n matrix

elems.U = unitaryfactory(n);

elems.e = euclideanfactory(n, 1);

manifold = productmanifold(elems);

problem.M = manifold;

problem.cost = @(X_struct) mycost(X_struct, M);

% solver choice is unclear; absent smoothness, would try pso or other

% derivative free solver, but it's likely better to entertain other cost

% functions to be used as proxies. This approximates the gradient (where

% possible) with finite differences, which is only good for quick prototyping

% in low dimensions. There are many ways to scale up from here.

% The solver can also take options for maximum iterations / time / etc.

X_struct = conjugategradient(problem);

U = X_struct.U;

lambda = 1.5 + .5*sin(X_struct.e);

X_matrix = U*diag(lambda)*U' % this is the X matrix found by the solver

function f = mycost(X_struct, M)

U = X_struct.U;

lambda = 1.5 + .5*sin(X_struct.e);

f = norm(diag(lambda.^(-1)) * (U'*M*U) * diag(lambda));

end

I hope this helps -- depending on interest, we could iterate here. Again, this problem does not naturally fit into the manopt framework for lack of smoothness, so it's likely that there are better tools for this problem out there.

Best,

Nicolas