Quadratic Objective with two Quadratic Constraints

Hii,

My programming knowledge is limited. I would like to solve a quadratic objective function with two quadratic constraints. I read and understood the concepts and examples given in https://in.mathworks.com/help/optim/ug/linear-or-quadratic-problem-with-quadratic-constraints.html?s_tid=mwa_osa_a, but I am having difficulty in implementing this for my problem.

Matrices A, B and C be symmetric positive semi-definite matrices. How to find the subject to the constraints and ? The documentation deals with a single quadratic constraint only. But my problem has two quadratic constraints. How can I code the second constraint as per the documentation in the link above?

Thank You..

[Matt J]

It seems like it would be a very difficult problem to solve for B~=C. Wth n=2, for example, the constrained region corresponds to the intersection of 2 ellipses, which is a discrete, non-connected set.

[Krishnendu K]

Is there any condition can I check whether the constrained region have an intersection?

I have generated the matrices A, B and C and attached the .mat files. Is there any possible way of optimisation?

[Matt J]

Even if they have an intersection, I don't see how that solves the problem. Do you have an initial guess close enough to the global solution that disconnected solution regions will be avoided?

Answer 1

[Torsten]

What about

function [y,yeq,grady,gradyeq] = quadconstr(x,B,C)
  y = [];
  yeq(1) = x.'*B*x - 1;
  yeq(2) = x.'*C*x - 1;
  if nargout > 2
    grady = [];
    gradyeq(:,1) = 2*B*x;    % Assumes B is symmetric, otherwise (B+B.')*x
    gradyeq(:,2) = 2*C*x;    % Assumes C is symmetric, otherwise (C+C.')*x
  end
end

[Krishnendu K]

Hii,

Thank You for the answer. How to create the hessian function for the same?

[Torsten]

My guess is

function hess = quadhess(x,lambda,A,B,C)
  hess = -2*A + lambda.eqnonlin(1)*2*B + lambda.eqnonlin(2)*2*C;
end

[Krishnendu K]

fun = @(x)quadobj(x,A);

nonlconstr = @(x)quadconstr(x,B,C);

x0 = zeros(N,1); % Column vector

options = optimoptions(@fmincon,'Algorithm','interior-point',...

'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true);

[x,fval,eflag,output,lambda] = fmincon(fun,x0,...

[],[],[],[],[],[],nonlconstr,options);

In this code I tried to initialize x with zeros, but got the resulting optimized vector as zero vector only. Is there any way I can initialize x?

[Torsten]

I hope you defined

function [y,grady] = quadobj(x,A)
  y = -x.'*A*x ;
  if nargout > 1
    grady = -2*A*x ;
  end
end

If you want us to check for possible problems, we need A, B and C.

[Krishnendu K]

Still there is problem with the initial values of x. Can you pls check whether the code works fine for the attached A, B and C matrices of n=6 size?

Amat = load("Matrix_A.mat");
A = Amat.P1

A = 
   0.1667 + 0.0000i   0.0625 + 0.0361i  -0.0208 - 0.0361i  -0.0000 - 0.0000i  -0.0208 + 0.0361i   0.0625 - 0.0361i
   0.0625 - 0.0361i   0.9167 + 0.0000i   0.2500 + 0.1443i  -0.2708 - 0.0361i   0.0000 - 0.0000i  -0.0833 + 0.1443i
  -0.0208 + 0.0361i   0.2500 - 0.1443i   0.9167 + 0.0000i   0.0625 + 0.0361i  -0.3333 - 0.1443i   0.0000 - 0.0000i
  -0.0000 + 0.0000i  -0.2708 + 0.0361i   0.0625 - 0.0361i   0.6667 + 0.0000i   0.0625 + 0.0361i  -0.0208 - 0.0361i
  -0.0208 - 0.0361i   0.0000 + 0.0000i  -0.3333 + 0.1443i   0.0625 - 0.0361i   0.9167 + 0.0000i   0.2500 + 0.1443i
   0.0625 + 0.0361i  -0.0833 - 0.1443i   0.0000 + 0.0000i  -0.0208 + 0.0361i   0.2500 - 0.1443i   0.4167 + 0.0000i

Bmat = load("Matrix_B.mat");
B = Bmat.Lambda_n

B = 6×6

1.0000 0 -0.5000 0 0 0 0 1.0000 0 -0.5000 0 0 -0.5000 0 1.0000 0 -0.5000 0 0 -0.5000 0 1.0000 0 -0.5000 0 0 -0.5000 0 1.0000 0 0 0 0 -0.5000 0 1.0000

Cmat = load("Matrix_C.mat");
C = Cmat.Gamma_n

C = 
   0.3333 + 0.0000i  -0.1250 - 0.0722i  -0.0417 - 0.0722i   0.0000 + 0.0000i  -0.0417 + 0.0722i  -0.1250 + 0.0722i
  -0.1250 + 0.0722i   1.3333 + 0.0000i   0.2500 + 0.1443i  -0.0417 - 0.0722i   0.0000 + 0.0000i   0.0833 - 0.1443i
  -0.0417 + 0.0722i   0.2500 - 0.1443i   1.3333 + 0.0000i  -0.1250 - 0.0722i   0.0833 + 0.1443i   0.0000 + 0.0000i
   0.0000 - 0.0000i  -0.0417 + 0.0722i  -0.1250 + 0.0722i   0.3333 + 0.0000i  -0.1250 - 0.0722i  -0.0417 - 0.0722i
  -0.0417 - 0.0722i   0.0000 - 0.0000i   0.0833 - 0.1443i  -0.1250 + 0.0722i   1.3333 + 0.0000i   0.2500 + 0.1443i
  -0.1250 - 0.0722i   0.0833 + 0.1443i   0.0000 - 0.0000i  -0.0417 + 0.0722i   0.2500 - 0.1443i   1.3333 + 0.0000i

N =size(A,1);
fun = @(x)quadobj(x,A);
nonlconstr = @(x)quadconstr(x,B,C);
x0 = zeros(N,1); % Column vector
options = optimoptions(@fmincon,'Algorithm','interior-point',...
    'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true);
[x,fval,eflag,output,lambda] = fmincon(fun,x0,...
    [],[],[],[],[],[],nonlconstr,options);

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance but constraints are not
satisfied to within the value of the constraint tolerance.


Consider enabling the interior point method feasibility mode.

function [y,grady] = quadobj(x,A)
  y = -x.'*A*x ;
  if nargout > 1
    grady = -2*A*x ;
  end
end
function [y,yeq,grady,gradyeq] = quadconstr(x,B,C)
  y = [];
  yeq(1) = x.'*B*x - 1;
  yeq(2) = x.'*C*x - 1;
  if nargout > 2
    grady = [];
    gradyeq(:,1) = 2*B*x;    % Assumes B is symmetric, otherwise (B+B.')*x
    gradyeq(:,2) = 2*C*x;    % Assumes C is symmetric, otherwise (C+C.')*x
  end
end

[Torsten]

Matrices A and C have complex entries. Something like max(x.'*A*x) does not exist for complex numbers.

[Krishnendu K]

In my case its not tranpose, its hermitian transpose, so i have used x'*A*x as the quadratic form. Thank you for your time and effort. Let me check where I have made the mistake in code. I am not sure whether this can be solved by using optimisation.

[Torsten]

x is not the problem, A and C are the problem.

I am not sure whether this can be solved by using optimisation.

What do you mean by "this" ?

[Krishnendu K]

What I meant is I tried to solve it by using eigen value decomposition. But the area of optimisation is entirely new to me. When I apply these constraints the problem cannot be solved by eigen value decomposition method and I have to rely entirely on convex optimisation which is new to me..

[Torsten]

What I meant is I tried to solve it by using eigen value decomposition.

You formulated an optimization problem. Obviously, the problem is not well-posed because the objective function is not real-valued.

What is the underlying problem that you tried to solved via this optimization formulation ? (You again talk of "I tried to solve it", but you don't explain what "it" is).