Matrix Eigenvalue/Positive definiteness optimization problem

Question

I have encountered a matrix optimization problem from a system of ODE's, however I am not very familiar with optimization.

Let $A(\lambda)$ be a $4 \times 4$ matrix dependent on a value $\lambda >0$ (some entries are constants, others are scalar multiples of $\lambda$) where my goal is to maximize this value of $\lambda$. However, there are some additional constraints:

1. the symmetric part of matrix the matrix product $-CA(\lambda)$ is positive definite (or equiv. strictly positive eigenvalues). Where the symmetric part is defined as $M_{symmetric} = \frac{1}{2}(M+M^{T})$

2. $C$ symmetric

3. $C$ preserves sign (Non-negative)

Here $C$ is a $4 \times 4$ matrix, which is to be chosen freely as long as it satisfies these constraints (the goal is to find a smart choice of $C$, maximizing this value of $\lambda$).

I have made this work for $A$ a $2 \times 2$ matrix in Mathematica (only 3 variables in $C$ using NMaximize[]). However for large $C$ this doesn't work anymore. Is there any software to tackle this problem with? Do you know similar problems with known steps to take? Any guidance in direction would be very welcome!

Answer 1

This seems to be a Bilinear Matrix Inequality (BMI) problem, which is non-convex and difficult. It is a BMI, rather than an easier LMI (Linear Matrix Inequality), because of the product of variables $C$ and $\lambda$. It can be formulated and attempted to be solved with a BMI-capable solver, such as PENLAB (free) or PENBMI (better, but not free), which can be used under YALMUP (under MATLAB). The formulation is not so difficult, but the solution might be.

The YALMIP code is something like:

C = sdpvar(4,4); % symmetric matric variable by default
lambda = sdpvar; % scalar variable
A = ...  supply the code to build up A from lambda
Constraints = [lambda >= 0, C(:) >= 0, -C*A-A'*C - min_eig*eye(4) >= 0];
optimize(Constraints,-lambda,sdpsettings('solver','penlab'))

where min_eig is chosen as a small positive number (minimum eigenvalue of -C*A-A'*C to ensure strict positive definiteness).

Given that this is a difficult nonlinear non-convex problem, supplying a non-default starting (initial) value for C and lambda might help a lot. That can be done with YALMIP's assign command and 'usex0',1 in sdpsettings. At best you will get a local optimum, which may not be globally optimal, The local optimum obtained might depend on the starting value.

Edit: Actually you could attempt to get a globally optimal solution, but may not succeed, using YALMUIP's BMIBNB branch and bound global optimization solver, using PENLAB (or better yet, PENBMI) as uppersolver (i.e. to solve local optimization problems spawned by BMIBNB). Change the optimize command to

optimize(Constraints,lambda,sdpsettings('solver','bmibnb','bmibnb.uppersolver','penlab'))

or better yet, specify penbmi rather penlab, if penbmi is available.

Even if BMIBNB does not succeed in finding the globally optimal solution (within some tolerance), it may find a lcally optimal solution as good, or possibly better, than what would have been found if only using penlab 9or penbmi) as a local solver.