A is a sp matrix, M is a PSD pp matrix, and I want to minimize det(AM^-1A^T). I believe AM^-1A^T can be formulated as a LMI, but the minimization of det (or equivalently of det_rootn) of R>=AM^-1A^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_endWhat 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+b11x1+b12x2;
y2 = a2+b21x1+b22x2;
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 = nmodnl^ndim;
np = nmod+nmodndim;
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)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);
%
% 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+V2lambda(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} …
diag([lambda(1)ones(1,np)/s1, lambda(2)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)));
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.