How to add the additional constraint on solution (and the optimization/constraint variables) that they need to be contained in the span/linear combination of a group of vectors?
I’m especially interested in the case where i have a hierarchy of constraints that are connected.
You will have to be more specific. Perhaps what you want can be done with affine (linear) constraints?
As for connected hierarchy of connected constraints, I have no idea what you’re talking about. If they are not convex, your only hope is if you can use MIDCP capability.
More precise explanation:
I define my optimization variables as following:
variable S(n,n,N+1) semidefinite;
minimize(trace(S(:,:,N+1)));
subject to
-S(:,:,1) <=PT_rho;
PT_rho<= S(:,:,1);
for i = 2:N+1
-S(:,:,i) <=PartialTranspose(S(:,:,i-1));
PartialTranspose(S(:,:,i-1))<=S(:,:,i);
end
I want to modify this code by adding a constraint, so that I force my S’s to be diagonalizable w.r.t. a certain choice of basis.
Meaning I can write` S=\sum_{i,j}^{d} x_{i,j} *(A_{j}\otimes B_{j}). Where A,B are some collection of matrices.
Is that possible with CVX, or do I need to use another package?
If the x_{i,.j} are variables, I think that is just an affine (linear) constraint.
So declare
variable x(n,n)
then add the constraint.
Does that do what you want?
But here variable x(n,n) is a matrix, so does it need to be defined properly. The idea is that i dont know them apriori, I only know what the different combinations of A_{i} and B_{j} and then the question becomes if there exists a sequence of constants such that i can rewrite them.
And have would one write such a sum within the cvx environment, does the subject to allow for nested for-loops?
I’m quite new to cvx, and am learning it myself, that is why I’m asking!
it’s your problem, not mine, so I don’;t really understand what you want. if you only need the existence of a matrix of x such that the constraint n S is satisfied, you should be able to do it. if you mean something else, i don’t know what that is.
CVX does allow the use of for loops, including nested for loops, provided that the indices in the for loops do not involve CVX variables or expressions.
You can build up CVX expressions using for loops. However, vectorization which eliminates or minimizes the number of levels of for loops allows the CVX modeling (processing) to be faster. There are many examples in this forum of vectorization and of one or more levels of nested for loops.
Hi, I think i modified the constraint such that it can be written as a linear combination of certain matrices(they are created in another function) My problem now is that the answer i now get is smaller then if i dont have the subspace constraint. I dont know if I have enforced the subspace constraint as i should do in CVX. (I’ve also tried using “==” ass well as some other variations of the code below )
Code with the subspace constraint:
function [opt, hierarchy, mu] = l_hierarchy_subpace_constraint(N, rho)
PT_rho = PartialTranspose(rho);
d=size(rho,1);
% Generate P matrices
P_array = generate_P_array();
% Start CVX optimization
cvx_begin quiet
% Define variables
variable a(N+1,d);
expression S(d,d,N+1);
minimize(a(N+1,1));
subject to
P_value=zeros(64,64);
for n=1:64
P_value=P_value+((1/(2^3))*a(1,n)*P_array(:,:,n));
end
S(:,:,1)=P_value;
% Enforce -S_1 <= PT_rho <= S_1
-S(:,:,1) <= PT_rho;
PT_rho <= S(:,:,1);
for j=2:N+1
P_value=zeros(64,64);
for n=1:64
P_value=P_value+((1/(2^3))*a(j,n)*P_array(:,:,n));
end
S(:,:,j)=P_value;
-S(:,:,j)<= PartialTranspose(S(:,:,j-1));
PartialTranspose(S(:,:,j-1)) <= S(:,:,j);
end
cvx_end
% Output results
opt = cvx_optval;
hierarchy= S;
mu = a;
end
First of all, don’t use quiet until you are sure your program is working properly. That way you can see all the solver and CVX output.
Are you saying that you have two programs, which are identical except for the addition of some (subspace constraints) in one of them; and that the program with the added constraints has a better (smaller) objective value than the other, and by more than just solver tolerance? If so, something must be amiss with your description or the execution of your program; and careful examination of solver and CVX output is called for. And if so, whether or not your constraints are “correct” (relative to what you want), something would be wrong. Caveat; if you are calling this in some iterative way such as SCA, then it would be possible that adding constraints could result in a better objective value, because one or both programs might not be finding the global minimum.
= is assignment for an expression.
== is an equality constraint.
There are no equality constraints in the code you posted. The only constraints are inequality constraints.
It is not apparent to me from your description where you tried using ==.
Your posts continue to be so vague as to make it difficult to help you.
Preformatted textHi Mark, first of all i want to say thank you for helping me. I will try to be more clear/precise. I have the following cvx, code that works wonderfully:
function[opt,hierarchy]=l_hierarchy(N,rho)
n = size(rho,1);
PT_rho=PartialTranspose(rho);
cvx_begin sdp
variable S(n,n,N+1) semidefinite;% This is a 3D slice where S(:,:,i+1) is S_i
minimize(trace(S(:,:,N+1)));
subject to
-S(:,:,1) <=PT_rho;
PT_rho<= S(:,:,1);
for i = 2:N+1
-S(:,:,i) <=PartialTranspose(S(:,:,i-1));
PartialTranspose(S(:,:,i-1))<=S(:,:,i);
end
cvx_end
opt=cvx_optval;
hierarchy=S;
end
My goal is to add another constraint to my optimization problem. This being a subspace constraint. I am able to to this in smaller dimentions, and where i don’t save alle the constraints-matrices (S-matrices) from my original problem. I use this code here:
function [opt,S_N,S_n] = l_subspace_d2(N,rho)
I=eye(4,4);
sx=[0,1;1,0];
sy=[0,-1i;1i,0];
sz=[1,0;0,-1];
PT_rho=PartialTranspose(rho);
cvx_begin
% Define variables
variable a(N+1,4);
minimize(trace(a(N+1,1)*I+a(N+1,2)*kron(sx,sx.')+a(N+1,3)*kron(sy,sy.')+a(N+1,3)*kron(sz,sz.')));
subject to
% Enforce -S_1 <= PT_rho <= S_1
-(a(1,1)*I+a(1,2)*kron(sx,sx.')+a(1,3)*kron(sy,sy.')+a(1,4)*kron(sz,sz.')) <= PT_rho;
PT_rho <= a(1,1)*I+a(1,2)*kron(sx,sx.')+a(1,3)*kron(sy,sy.')+a(1,4)*kron(sz,sz.');
for i=2:N+1
-(a(i,1)*I+a(i,2)*kron(sx,sx.')+a(i,3)*kron(sy,sy.')+a(i,4)*kron(sz,sz.'))<=PartialTranspose(a(i-1,1)*I+a(i-1,2)*kron(sx,sx.')+a(i-1,3)*kron(sy,sy.')+a(i-1,4)*kron(sz,sz.'));
PartialTranspose(a(i-1,1)*I+a(i-1,2)*kron(sx,sx.')+a(i-1,3)*kron(sy,sy.')+a(i-1,4)*kron(sz,sz.'))<=a(i,1)*I+a(i,2)*kron(sx,sx.')+a(i,3)*kron(sy,sy.')+a(i,4)*kron(sz,sz.');
end
cvx_end
opt = cvx_optval;
S_N = a(N+1,1)*I+a(N+1,2)*kron(sx,sx.')+a(N+1,3)*kron(sy,sy.')+a(N+1,4)*kron(sz,sz.');
S_n = a(N,1)*I+a(N,2)*kron(sx,sx.')+a(N,3)*kron(sy,sy.')+a(N,4)*kron(sz,sz.');
end
My problem comes when I want to generalize to higher dimentional matrices, which i create with the following function:
function P_array = generate_P_array()
I = [1, 0; 0, 1];
sx = [0, 1; 1, 0];
sy = [0, -1i; 1i, 0];
sz = [1, 0; 0, -1];
pauli_set = {I, sx, sy, sz};
P_array = zeros(64,64, 64);
i = 1;
for n1 = 1:4
for n2 = 1:4
for n3 = 1:4
P = kron(kron(pauli_set{n1}, pauli_set{n2}), pauli_set{n3});
P_t = kron(kron(pauli_set{n1}.', pauli_set{n2}.'), pauli_set{n3}.');
P_array(:,:,i) = kron(P, P_t); % Tensor with transposed version
i = i + 1;
end
end
end
end
My attempt at a code that tries to generalize the optimisation problem:
function [opt, hierarchy, mu] = l_subpace_gen(N, rho)
PT_rho = PartialTranspose(rho);
% Generate the Pauli basis expansion
P_array = generate_P_array(); % P_array is d × d × d (64×64×64)
% Start CVX optimization
cvx_begin sdp
% Define optimization variables
variable a(N+1, 64); % Coefficients for the Pauli basis
expression S(64, 64, N+1); % Store the S matrices
% Construct S matrices in the Pauli basis
for j = 1:N+1
%X=reshape(a(j, :), [1, 1, 64]);
%S(:,:,j) = sum(P_array .* X, 3);
S(:,:,j) = sum(P_array .* reshape(a(j, :), [1, 1, 64]), 3);
end
% Objective function
minimize(trace(S(:,:,N+1)));
% Constraints enforcing subspace structure
S(:,:,1) >= -PT_rho;
S(:,:,1) <= PT_rho;
for j = 2:N+1
S(:,:,j) >= -PartialTranspose(S(:,:,j-1));
S(:,:,j) <= PartialTranspose(S(:,:,j-1));
end
cvx_end
% Output results
opt = cvx_optval;
hierarchy = S;
mu = a;
end
The main error i encounter is this:
Error using .* (line 46)
Matrix dimensions must agree.
Error in l_subpace_gen (line 18)
S(:,:,j) = sum(P_array .* X, 3);
Even though when i test the following:
P=generate_P_array();
a=rand(4,64);
S=zeros(64,64,4);
for j = 1:4
%X=reshape(a(j, :), [1, 1, 64]);
%S(:,:,j) = sum(P_array .* X, 3);
S(:,:,j) = sum(P .* reshape(a(j, :), [1, 1, 64]), 3);
end
I dont get any errors. Thus it may seem that this is a problem with cvx, does declaring a and S as variable and expression make MATLAB view them differently? I know that they are not seen as double/int since I cant use “pagemtimes” function, this would also not be possible for MATLAB versions older then 2020(but I have 2024b).
Hope this clarifies everything, and that you may be able to help my predicament.
You still haven’t shown “clean” code which when executed results in your error message. So I am nos sure what you ran which resulted in the error
Error using .* (line 46)
Matrix dimensions must agree.
Error in l_subpace_gen (line 18)
S(:,:,j) = sum(P_array .* X, 3);
Based on the error message, it appears that P_array and X do not have the same dimensions.
Note that CVX does not use implicit expansion. So If either or both of P_array and X are CVX expressions or variables, their dimensions need to be identical, because implicit expansion is not allowed. Therefore, you must use repmat, reshape or something to get dimensions identical on each side of .*
Sorry, I see i have pasted the wrong error message, i run the following code
N=5;
[l_2,Z]=l_subpace_gen(N,gamma);
l_2
where gamma is a prechoosen 64x64 matrix. And l_supace_gen() is the function defined above. I get the following error:
Error using .* (line 46)
Matrix dimensions must agree.
Error in l_subpace_gen (line 17)
S(:,:,j) = sum(P_array .* reshape(a(j, :), [1, 1, 64]), 3);
Is this correct way of expressing clean code?
if P_array is 8 by 8 by 8 as your comment says, obviously it is incompatible with 1 by 1 by 64. You need to have identical dimensions for both. You can’t rely on implicit expansion. You can use repmat.
Would not P_array be a 6464 matrix since you get 2^{6} from the generate_P_array function, since you tensor 6 22 matrices?
It’s your code and your comments.
You can use whos on MATLAB variables to see their dimensions; and type CVX expressions at the command line (or in a function) to see their dimensions.
Running
disp(size(reshape(a(j, :), [1, 1, 64])))
disp(size(P_array))
yields the output
1 1 64
64 64 64
So the dimensions match in the 3rd dimension. Which i think is the one I multiply element wise and then sum over when running sum(P_array .* reshape(a(j, :), [1, 1, 64]), 3) or am I wrong here?
Yes, you are wrong. .* occurs before sum. Both sides of .*` must have identical dimensions.
But then why does running
P_array=generate_P_array();
a=rand(4,64);
S=zeros(64,64,4);
for j = 1:4
S(:,:,j) = sum(P_array .* reshape(a(j, :), [1, 1, 64]), 3);
end
in a non cvx environment, not result in the same error message?
Here i do get a 3d array S where S(:,:,j) are 64*64 matrices where i have multiplied each element of P_array with an element of a(j,:).
As I wrote, CVX does not allow implicit expansion.
https://blogs.mathworks.com/loren/2016/11/10/more_thoughts_about_implicit_expansion/ does not apply to CVX expressions or variables.
Thank you for all the help Mark!