Optimizing matrix inequalities over trace

I need to solve the following problem.

Given $n \times n$ Hermitian matrices $A\geq 0$ and $B_1, ~ B_2$ (need not be positive semidefinite), with $Tr(AB_1)<0~,Tr(AB_2)<0$ construct a Hermitian matrix $B$ such that $B\geq B_1$ and $B\geq B_2$, for which $\min Tr(AB)$ is achieved.

For Hermitian matrices $X$ and $Y$, we write $X\geq Y$ if and only if $X-Y\geq 0$ (i.e. we define a partial order). There is no partial order relations between $B_1$ and $B_2$, otherwise the solution would have been easier.

Using some earlier programmes, I have constructed $A$, $B_1,~B_2$ (as well as $B_3,~B_4,...$, finite number of such $B$'s). Needless to say, there are all huge matrices and was time consuming to construct. I really do not want to extract the matrices as data and use it in some other software, unless it is absolutely necessary. Can someone please help me to pointing out some method for solving in Mathematica itself?

If I have not tagged it correctly,I request the moderators to help. Feel free to ask me any question if you need any clarification.

1 Answer

Here is one approach. I have not tested it so there might be unforseen hitches. I'll refer to your matrices as matrixA, matrixB1, and matrixB2 respectively.

Define

minEig[mat1:{{_Real..}..},mat2:{{_Real..}..}] :=
Min[Eigenvalues[mat1-mat2]]
newB = Array[b, Dimensions[a]];
obj = Trace[matrixA.newB];
FindMinimum[{obj, minEig[newB,matrixB1]>=0, minEig[newB,matrixB2]>=0},
Flatten[newB]]


I'd be curious to find out if this actually works...

--- edit ---

Okay, I finally got back to this one. First thing to note is that, like every other Matheamtica programmer on the planet doing linear algebra, I managed to mix Trace for Tr. Okay, on to the chase.

Here is the basic input for the example provided in a comment to the original post.

matA = {{1, I}, {-I, 2}};
matB1 = {{1/2, 1/2 (I + Sqrt[3])}, {1/2 (-I + Sqrt[3]), -1}};
matB2 = {{1, 1/2 (-I - Sqrt[3])}, {1/2 (I - Sqrt[3]), -1}};


We set up our variables and result matrix.

n = Length[matA];
matX = Array[x, {n, n}];
matY = Array[y, {n, n}];
Do[matX[[i, j]] = matX[[j, i]]; matY[[i, j]] = -matY[[j, i]];
, {i, 2, n}, {j, 1, i - 1}];
Do[matY[[i, i]] = 0, {i, n}];
matB = matX + I*matY;
vars = Variables[{matX, matY}];


Now we need an objective function and something to handle the (matrix) nonnegativity constraints.

obj = Expand[Tr[matA.matB]];
minEig[mat1 : {{_?NumberQ ..} ..}, mat2 : {{_?NumberQ ..} ..}] :=
Min[Eigenvalues[mat1 - mat2]]


We're off to the races.

FindMinimum[{obj, minEig[matB, matB1] >= 0,
minEig[matB, matB2] >= 0}, vars]

(* Out[199]= {0.886000815253, {x[1, 1] -> 2.53487711894,
x[1, 2] -> 0.30408162278, x[2, 2] -> -0.0636708567376,
y[1, 2] -> -0.760767295104}} *)


I do not know how good a result this is, but at least it might provide a start in terms of coding something useful for the task at hand.

--- end edit ---

• Thank you for your help. I am trying to implement it with simple examples. I guess we can replace Real by Complex. Isn't it? My matrices are not real matrices, they are Hermitian. I am working on the program. I shall keep on informing you as it progresses. I am also concerned that FindMinimum may take a long time to run and give output. – RSG Mar 18 '15 at 16:03
• Yes you can use Complex but you may need to extract Re on eigenvalues and obj in order to handle numerical artifact complex parts. – Daniel Lichtblau Mar 18 '15 at 16:15
• There is error. I guess it is due the following things. The matrix newB is not defined as Hermitian matrix. Hence, we are getting all types of complex numbers as outputs of minEig and Tr. Is it possible to restrict our search only on the Hermitian matrices? – RSG Mar 18 '15 at 16:27
• Yeah, I forgot about that. Not too hard to define as Hermitian but it will mean splitting into explicitly real and imaginary parts. So you'd need newBRe, newBIm, run a double loop that sets lower equal upper part on the first, and negative of upper on the second. Then have newB=newBRe+newBIm. Also the variable list will require something like Variables[newB] so as to avoid duplicates. – Daniel Lichtblau Mar 18 '15 at 19:46
• By the way, do you have a smallish specific example to work with? – Daniel Lichtblau Mar 18 '15 at 20:44