minimum value optimization of matrix with constraints

i have a matrix A of size n×n such that all of its entries are known excpt m number of entries are unknown.

then how i will calculate the minimum of follwing objective function subject to A are symmetric positive semidefinite (i.e. all eigenvalue of A are greater than or equal to 0)

obj = @(A) 0.5 * norm(A, 'fro')^2
constraint = @(A) (all (eig(X_new>=0)))

Answer 1

[Bruno Luong]

Straighforward implementation, it seems more reliable than Matt's square-root parametrization.

The square-root makes the gradient close to 0 close to the boundary of the admissible set, and it can create somewhat difficulty of the convergence by gradient descend methods such as fmincon.

while true
    n=3;
    [Q,~]=qr(randn(n));
    d=[rand(n-1,1);0];
    A=Q*diag(d)*Q';
    if any(A < 0, 'all')
        break
    end
end
A

A = 3×3

0.3736 -0.1307 -0.3455 -0.1307 0.5302 0.1397 -0.3455 0.1397 0.3202

Iunknown=A<0;
Adecimation = A;
Adecimation(Iunknown) = 0

Adecimation = 3×3

0.3736 0 0 0 0.5302 0.1397 0 0.1397 0.3202

i = 1:n;
Iunknown = Iunknown & i>=i'; % work only unknown in upper triangular part
x0=zeros(nnz(Iunknown),1);
norm(A,'fro')

ans = 0.9139

x=fmincon(@(x) fronorm(x,Iunknown,A),x0,[],[],[],[],[],[],@(x) nonlcon(x,Iunknown,A));

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

Afmincon = CompletionA(x, Iunknown, Adecimation),

Afmincon = 3×3

0.3736 0.2020 -0.1010 0.2020 0.5302 0.1397 -0.1010 0.1397 0.3202

c=nonlcon(x,Iunknown,Adecimation);
% Check validity of the solution
OK = norm(Afmincon,'fro')-norm(A,'fro') < 1e-6 || ...
    c >= 0;
if ~OK
    fprintf('fails to find global minima solution\n')
    norm(Afmincon,'fro')
    norm(A,'fro')
end

function A = CompletionA(x, Iunknown, Adecimation)
A=Adecimation;
A(Iunknown) = x;
A = 0.5*(A+A.');
end

function f=fronorm(x,Iunknown,Adecimation);
A = CompletionA(x, Iunknown, Adecimation);
f = sum(A.^2,"all");
end

function [c,ceq]=nonlcon(x,Iunknown,Adecimation)
ceq=[];
A = CompletionA(x, Iunknown, Adecimation);
c=min(eig(A));
end

Answer 2

[Matt J]

Perhaps as follows:

%L=unknown parameter matrix (nxn) where A=L.'*L
%Iknown=indices of known elements  of A (mx1)
%Aknown=known A values (mx1)


L0=Init(Iknown,Aknown,n);

L=fmincon(@(L) 0.5*norm(L*L.', 'fro')^2,  L0,[],[],[],[],[],[],...
          @(L) nonlcon(L,Iknown,Aknown), options );

A=L*L.'; %optimal A

function [c,ceq]=nonlcon(L,Iknown,Aknown)

    c=[];
    
    A=L*L.';
    
    ceq=A(Iknown)-Aknown;

end

function L0=Init(Iknown,Aknown,n)

  A=zeros(n);
  A(Iknown)=Aknown;

  A=A-(min(eig(A))-0.01)*eye(n);

  L0=chol(A);

end

[Bruno Luong]

I like Matt's formulation

But I don't know how to initialize L0, it looks like fmincon fails to find global minima now and then.

while true
    n=3;
    [Q,~]=qr(randn(n));
    d=[rand(n-1,1);0];
    A=Q*diag(d)*Q';
    if any(A < 0, 'all')
        break
    end
end
A

A = 3×3

0.2827 -0.0446 -0.0496 -0.0446 0.0070 0.0064 -0.0496 0.0064 0.8655

Iknown=A>=0;
L0=zeros(n);
Aknown = A(Iknown);
norm(A,'fro')

ans = 0.9154

obj = @(L) 0.5*norm(L*L.', 'fro')^2;
L=fmincon(obj,L0,[],[],[],[],[],[],@(L) nonlcon(L,Iknown,Aknown) );

Warning: Matrix is singular to working precision.

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.

Afmincon=L*L.' %optimal A

Afmincon = 3×3

0.1331 0.0033 0.4075 0.0033 0.0001 0.0101 0.4075 0.0101 1.2474

OK = norm(Afmincon,'fro')-norm(A,'fro') < 1e-6;
if ~OK
    fprintf('fails to find global minima solution\n')
    norm(Afmincon,'fro')
    norm(A,'fro')
end

fails to find global minima solution

ans = 1.3806

ans = 0.9154

function [c,ceq]=nonlcon(L,Iknown,Aknown)
%L=unknown variables (nxn)
%Iknown=indices of known matrix elements (mx1)
%Aknown=known matrix values 9mx1)
c=[];

A=L*L.';

ceq=A(Iknown)-Aknown;
end

[Matt J]

This initialization scheme seems to work well:

k=0;
N=300;
while k<N

    OK=doOptim();
    
    if OK
     k=k+1;
    else
        break;
    end

end
disp("Success rate = " + 100*k/N + "%")

Success rate = 100%

function OK=doOptim()

    while true
        n=3;
        [Q,~]=qr(randn(n));
        d=[rand(n-1,1);0];
        A=Q*diag(d)*Q';
        if any(A < 0, 'all')
            break
        end
    end
    
    Iknown=A>=0;
    
    Aknown = A(Iknown);
    
    L0=Init(Iknown,Aknown,n);

    obj = @(L) 0.5*norm(L*L.', 'fro')^2;
    opts=optimoptions('fmincon','Display','none');
    warning off
    L=fmincon(obj,L0,[],[],[],[],[],[],@(L) nonlcon(L,Iknown,Aknown) ,opts);
    warning on
    Afmincon=L*L.'; %optimal A

    OK = norm(Afmincon,'fro')-norm(A,'fro') < 1e-6;

end

function [c,ceq]=nonlcon(L,Iknown,Aknown)
%L=unknown variables (nxn)
%Iknown=indices of known matrix elements (mx1)
%Aknown=known matrix values 9mx1)
c=[];

A=L*L.';

ceq=A(Iknown)-Aknown;
end

function L0=Init(Iknown,Aknown,n)

  A=zeros(n);
  A(Iknown)=Aknown;

  A=A-(min(eig(A))-0.01)*eye(n);

  L0=chol(A);

end

[Bruno Luong]

+1