A is a s*p matrix, M is a PSD p*p matrix, and I want to minimize det(A*M^-1*A^T). I believe A*M^-1*A^T can be formulated as a LMI, but the minimization of det (or equivalently of det_rootn) of R>=A*M^-1*A^T is not possible. In practice this has some similarity with the function matrix_frac except here x is a matrix and we are interested in a function of the resulting matrix.

Any ideas or suggestions are welcome.

If s = 1,` minimize(matrix_frac(A',M))`

If s = p, `minimize(det_inv(M))`

(you can leave out the factor of `det(A)^2`

)

Otherwise, `?`

Have you proven it is convex? You can see the code for `matrix_frac`

and `det_inv`

as well as `det_rootn`

in `functions/@cvx`

. I don’t see any obvious modifications where you could use `A`

or `A'`

in the off-diagonal terms in place of what’s there and end up with a relevant LMI.

Thank you, Mark. Yes, 1<s<p. Also, it is convex. I just ploted it. In simplest case it can be reduced to a single decision variable. Finally, about your suggestion, that was also the idea I was exploiting indeed. I send a piece of code below. However, I am getting an error which is expected … det_rootn is indeed a concave function and I am using in a minimization. Do you see a way to overcome it? I am probably doing a very naif error.

Disciplined convex programming error:

Cannot minimize a(n) concave expression.

```
cvx_precision high
variable lambda(nm)
variable Rmat(nr,nr) symmetric
Vmat = V1*lambda(1)/s1+V2*lambda(2)/s2;
minimize (det_rootn(Rmat))
subject to
[Vmat, transpose(Amat); Amat, Rmat] == semidefinite(npar+nr);
sum(lambda) == 1;
lambda >= 0;
cvx_end
```

What is your proof of convexity? Is it constructive (in DCP sense)?

I just varied the decision variable in its domain which is [0,1], and for each value computed the function. Then, plotted the function and it is a typical convex function. This function corresponds to D_A-optimal design of experiments and is known to be convex, see Atkinson, Donev and Tobias (2007).

Huh? `M`

is a p by p matrix, having `p*(p+1)/2`

unique elements. The function needs to be jointly convex in all matrix elements. it is not one dimensional. This is the first you mentioned of a [0,1] domain.

I have not, and still have not, seen that reference.

Sorry for only now mentioning the domain. I am interested in a formulation where the decision variables are positive. In this case and for the purpose of simplyfing we can reduce the number of decision variables to 1 and its domain to [0,1]. Matrix M is a Fisher Information Matrix (SDP). In this particular case M is a diagonal matrix, but generally it can be a full matrix having p*(p+1)/2 unique elements, as you say.

In case you want I can send you the code I am using for testing.

Now you are telling us the rather significant seeming fact that M is diagonal. Y

You say that you can reduce the number of decision variables to one. Does that mean that in fact `M`

is a scalar multiple of the identity matrix? Why don’t you show us your one variable formulation?

Sure. I can share the code for this example. Thanks for your help.

function example1

%

% nmod - number of models

% ndim - number of factors x

% nc - number of total candidate points for nmod models

% np - number of parameters

% i_model - number of candidate points per model

clear all;

global Amat

%

syms t w k x1 x2 a1 a2 b11 b12 b21 b22

y1 = a1+b11*x1+b12*x2;

y2 = a2+b21*x1+b22*x2;

h11 = simplify([diff(y1,a1); diff(y1,a2); diff(y1,b11); diff(y1,b12); …

diff(y1,b21); diff(y1,b22)]);

h12 = simplify([diff(y2,a1); diff(y2,a2); diff(y2,b11); diff(y2,b12); …

diff(y2,b21); diff(y2,b22)]);

%

nmod = 2;

ndim = 2;

nl = 2;

nc = nmod*nl^ndim;
np = nmod+nmod*ndim;

i_model = nc/nmod;

%

%

%

% initial grid

low_x(1:ndim) = [-1, -1];

up_x(1:ndim) = [1, 1];

dx(1:ndim) = [2, 2];

%

% A-matrix

Amat = double([1, -1, 0, -0, 0, 0; …

0, 0, 1, -1, 0, 0]);

%Amat = [1, -1, 1, -1, 0, 0];

%

%

for i=1:ndim

nx(i) = (up_x(i)-low_x(i))/dx(i)+1;

end

nxmax = max(nx);

ncomb = prod(nx);

%

% construct candidate points

[x,nz] = allcandidates(low_x,up_x,dx,nx);

%

dh = zeros(np,nc);

fim_ind = zeros(nmod,ncomb,np,np);

for m=1:nmod

for i=1:ncomb

ir1 = x(i,1);

jr1 = x(i,2);

mr1 = (m-1)

*ncomb+i;*

if m==1

dh(1:np,mr1) = subs(h11,{x1,x2},{ir1, jr1});

else

dh(1:np,mr1) = subs(h12,{x1,x2},{ir1, jr1});

end

fim_ind(m,i,1:np,1:np) = double(dh(1:np,mr1)

%

% prepare to solve conic problem

%

% initialization of the solver

cvx_setup

%

%

s1 = 1.0;

s2 = 0.5;

alls1 = [0.2:0.2:1];

alls2 = fliplr(alls1(1:end-1));

alls = [alls1, 1./alls2];

ns = size(alls,2);

savex = zeros(3,ns,nmod);

saveo = zeros(3,ns);

for i=1:1 % criteria (only D_A; i.e, ind_criterion == 1 in function optimal_design)

for j=1:ns

s2 = alls(j);

[xf, obj_dec, nw] = optim_design(V1, V2, nmod, ndim, np, nc, s1, s2, i);

saves(j,1) = s1;

saves(j,2) = s2;

savex(i,j,1:nmod) = xf(1:nmod,2);

saveo(i,j) = obj_dec;

end

end

%

% print results

%file1 = fopen(‘Ccovar2lev.txt’,‘w’);

%str = ‘\’;

%for j=1:ns

% fprintf(file1,’%6.4f & %6.4f && ‘, saves(j,1:2));

% for i=1:3

% fprintf(file1,’%6.4f & %6.4f ‘, savex(i,j,1));

% if i==3

% fprintf(file1,’%8.4f %s \n’, saveo(i,j), str);

% else

% fprintf(file1,’%8.4f && ', saveo(i,j));

% end

% end

%end

%fclose(file1);

end

%% solve SDP problems

%

function [xf, obj_dec, nw] = optim_design(V1, V2, nm, nx, npar, nq, s1, s2, ind_criterion)

%

% solve the SDP problem for a given criterion

global Amat

np = ceil(nq/nm);

nr = size(Amat,1);

addpath([cvx_where,‘functions’,filesep,‘vec_’]);

savepath

cvx_solver Mosek

%

if ind_criterion == 1 %% DA-optimality criterion

%

% solve the SDP problem for DA-optimality

cvx_begin sdp

cvx_precision high

variable lambda(nm)

variable Rmat(nr,nr) symmetric

Vmat = V1lambda(1)/s1+V2

if m==1

dh(1:np,mr1) = subs(h11,{x1,x2},{ir1, jr1});

else

dh(1:np,mr1) = subs(h12,{x1,x2},{ir1, jr1});

end

fim_ind(m,i,1:np,1:np) = double(dh(1:np,mr1)

*transpose(dh(1:np,mr1))/ …*

i_model);

end

end

fim_mod = zeros(np,np,nmod);

for m=1:nmod

fim_mod(1:np,1:np,m) = sum(fim_ind(m,1:ncomb,1:np,1:np),2);

end

V1 = fim_mod(1:np,1:np,1);

V2 = fim_mod(1:np,1:np,2);

%V1{1} = double(dh(1:np,1:i_model)/i_model);

%V2{1} = double(dh(1:np,i_model+1:2i_model)/i_model);i_model);

end

end

fim_mod = zeros(np,np,nmod);

for m=1:nmod

fim_mod(1:np,1:np,m) = sum(fim_ind(m,1:ncomb,1:np,1:np),2);

end

V1 = fim_mod(1:np,1:np,1);

V2 = fim_mod(1:np,1:np,2);

%V1{1} = double(dh(1:np,1:i_model)/i_model);

%V2{1} = double(dh(1:np,i_model+1:2

%

% prepare to solve conic problem

%

% initialization of the solver

cvx_setup

%

%

s1 = 1.0;

s2 = 0.5;

alls1 = [0.2:0.2:1];

alls2 = fliplr(alls1(1:end-1));

alls = [alls1, 1./alls2];

ns = size(alls,2);

savex = zeros(3,ns,nmod);

saveo = zeros(3,ns);

for i=1:1 % criteria (only D_A; i.e, ind_criterion == 1 in function optimal_design)

for j=1:ns

s2 = alls(j);

[xf, obj_dec, nw] = optim_design(V1, V2, nmod, ndim, np, nc, s1, s2, i);

saves(j,1) = s1;

saves(j,2) = s2;

savex(i,j,1:nmod) = xf(1:nmod,2);

saveo(i,j) = obj_dec;

end

end

%

% print results

%file1 = fopen(‘Ccovar2lev.txt’,‘w’);

%str = ‘\’;

%for j=1:ns

% fprintf(file1,’%6.4f & %6.4f && ‘, saves(j,1:2));

% for i=1:3

% fprintf(file1,’%6.4f & %6.4f ‘, savex(i,j,1));

% if i==3

% fprintf(file1,’%8.4f %s \n’, saveo(i,j), str);

% else

% fprintf(file1,’%8.4f && ', saveo(i,j));

% end

% end

%end

%fclose(file1);

end

%% solve SDP problems

%

function [xf, obj_dec, nw] = optim_design(V1, V2, nm, nx, npar, nq, s1, s2, ind_criterion)

%

% solve the SDP problem for a given criterion

global Amat

np = ceil(nq/nm);

nr = size(Amat,1);

addpath([cvx_where,‘functions’,filesep,‘vec_’]);

savepath

cvx_solver Mosek

%

if ind_criterion == 1 %% DA-optimality criterion

%

% solve the SDP problem for DA-optimality

cvx_begin sdp

cvx_precision high

variable lambda(nm)

variable Rmat(nr,nr) symmetric

Vmat = V1

*lambda(2)/s2;*

minimize (det_rootn(Rmat))

subject to

[Vmat, transpose(Amat); Amat, Rmat] == semidefinite(npar+nr);

sum(lambda) == 1;

lambda >= 0;

cvx_end

%

% optimum

obj_dec = det_rootn(Rmat);

elseif ind_criterion == 2

%

% solve the SDP problem for A-optimality

cvx_begin sdp

cvx_precision default

variable lambda(nm)

Vmat = Vall{1}…

minimize (det_rootn(Rmat))

subject to

[Vmat, transpose(Amat); Amat, Rmat] == semidefinite(npar+nr);

sum(lambda) == 1;

lambda >= 0;

cvx_end

%

% optimum

obj_dec = det_rootn(Rmat);

elseif ind_criterion == 2

%

% solve the SDP problem for A-optimality

cvx_begin sdp

cvx_precision default

variable lambda(nm)

Vmat = Vall{1}

diag([lambda(1)

*ones(1,np)/s1, lambda(2)*

Vall{1}’;

minimize (trace_inv(Vmat))

subject to

Vmat == semidefinite(npar);

sum(lambda) == 1;

lambda >= 0;

cvx_end

%

% optimum

obj_dec = trace_inv(Vmat);

elseif ind_criterion == 3

%

% solve the SDP problem for E-optimality

cvx_begin sdp

cvx_precision high

variables lambda(nm)

Vmat = Vall{1}…

*ones(1,np)/s2])*…Vall{1}’;

minimize (trace_inv(Vmat))

subject to

Vmat == semidefinite(npar);

sum(lambda) == 1;

lambda >= 0;

cvx_end

%

% optimum

obj_dec = trace_inv(Vmat);

elseif ind_criterion == 3

%

% solve the SDP problem for E-optimality

cvx_begin sdp

cvx_precision high

variables lambda(nm)

Vmat = Vall{1}

diag([lambda(1)*ones(1,np)/s1, lambda(2)

*ones(1,np)/s2])*…

Vall{1}’;

maximize (lambda_min(Vmat))

subject to

Vmat == semidefinite(npar);

sum(lambda) == 1;

lambda >= 0;

cvx_end

%

% optimum

obj_dec = lambda_min(Vmat);

end

%

% prunning

lambdaD = lambda;

[lambdaj] = find(lambda>-1e-4);

qt = sum(lambda(lambdaj),1);

nw = size(lambdaj,1);

x1z = [1:nm]’;

lambdaz1 = lambda(lambdaj)/qt;

xf = [x1z, lambdaz1];

return

end

%%

function [x,nz] = allcandidates(low_x,up_x,dx,nx)

%

% procedure that generates all candidate points

% O: x - set of candidates

% nz - #candidates

%

ndim = size(low_x,2);

nxmax = max(nx);

ncomb = prod(nx);

z = zeros(nxmax,ndim);

x = [];

x2 = [];

xnew = [];

for i=1:ndim

for j=1:nx(i)

z(j,i) = (j-1)

*dx(i)+low_x(i);*

end

end

x = allcomb(transpose(z(1:nx(1),1)), transpose(z(1:nx(2),2)));

nz = size(x,1);

k = 2;

for j=k+1:ndim

for i=1:k

if j==k+1

x2(1:nznx(3),1:2) = allcomb(x(1:nz,i),transpose(z(1:nx(j),j)));

end

end

x = allcomb(transpose(z(1:nx(1),1)), transpose(z(1:nx(2),2)));

nz = size(x,1);

k = 2;

for j=k+1:ndim

for i=1:k

if j==k+1

x2(1:nz

nznew = size(x2,1);

xnew(1:nznew,i) = x2(1:nznew,1);

xnew(1:nznew,j) = x2(1:nznew,2);

end

end

x = xnew;

nz = nznew;

k = k+1;

end

end

I have no idea what all that is or means, how it constitutes a single scalar variable formulation, or shows its convexity, and don’t have the energy to try to figure it out.

But I will offer the advice that you should not use `cvx_precision high`

with Mosek, or any other solver. Just stick with the default.

Thank you. Can I use cvx for problems of `max det(R)^{-1/k}`

; R is sdp matrix and k is its size? This problem is equivalent to that I have.

No.

`det(R)^{-1/k}`

is convex, and can be formulated as `inv_pos(det_rootn(R))`

So it can be minimized, not maximized, in CVX.