Minimizing the nuclear norm of a matrix

[Bhagyashri Telsang]

Hello,

The objective function I have is nuclear norm of a matrix X, along with one other norm term. I need to minimize this objective function and there are no constraints. The problem I am facing is that svd cannot be applied on sdpvar objects. So in order to go about this problem, I tried to supply a numerical initial value to X, using assign. When I do this, the value that I get after the optimization problem is solved is same as the value I had initialized it with. The code I am working on is as follows.

A_est = sdpvar(n,n);

B_est = sdpvar(n,1);

C_est = sdpvar(1,n);

D_est = sdpvar(1,1);

T_u_s = sdpvar(s,s);

T_y_s = sdpvar(s,s);

for i=1:s

    for j=1:s

        if i==j

            T_u_s(i,j) = D_est;

        elseif i<j

            T_u_s(i,j) = 0;

        elseif i==2 & j==1

            T_u_s(i,j) = C_est *  B_est;

        else

            T_u_s(i,j) = C_est * A_est^(i-2) * B_est;

        end

    end

end

   

for i=1:s

    for j=1:s

        if i==j

            T_y_s(i,j) = 0;

        elseif i<j

            T_y_s(i,j) = 0;

        else

            T_y_s(i,j) = C_est * A_est^(i-2) * B_est;

        end

    end

end

y_est = sdpvar(N,1);

Y_est = sdpvar(s,(N-s+1));

for i = 1:s

    for j = 1:N-s+1

        Y_est(i,j) = y_est(j+i-1);

    end

end

summation_of_norm = 0;

for i = 1:N

    summation_of_norm = summation_of_norm + (y(i)-y_est(i))^2;

end

summation_of_norm = summation_of_norm/N;

A0 = 0.73 * eye(n,n);

B0 = ones(n,1);

C0 = randn(1,n);

D0 = 0.2;

y_est0 = zeros(N,1);

assign(y_est, y_est0)

assign(A_est, A0)

assign(B_est, B0)               %ERROR: y_est is being approximated well but A,B,C,D remain the way they were initialised. 

assign(C_est, C0)               

assign(D_est, D0)

X = Y_est - T_y_s * Y;

constraint = [];

objective = norm_nuc(value(X))+ summation_of_norm ;    

DIAGNOSTIC = solvesdp(constraint,objective) ;

    

The matrix X is in terms of A,B,C,D which are optimization variables. Also, y_est is an optimization variable which is also initialized using the command 'assign'. y_est is being approximated well because it appears in the term summation_of_norm. But other optimization variables that appear only in X, do not change; they stay as initialized. 

I tried using 'usex0' but I get the error that the solver does not support warm-start through YALMIP and I have not been able to find a solver that supports 'usex0'. 

The main problem is minimizing the nuclear norm of X. It would be great if you could help me out with this. Thank you in advance! 

[Johan Löfberg]

First, note that all square matrices you have defined are symmetric. Probably not what you want.
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Tutorials.Basics

Nuclear norms are treated in YALMIP using norm(X,'nuclear') in a convex semidefinite programming context

What is the function norm_nuc?

You say that you have solved the problem but doubt the answer. Did you actually look at the diagnostic from solvesdp? Did a solver actually run?

However, you are setting up an absolutely horrible nonlinear nonconvex model. Trying to use YALMIP to solve this problem is the wrong tool. To be honest, I think you have to take a step back and try to come up with better model altogether. Really, you have a nasty nonconvex problem, so why do you bother with the introduction of the nuclear norm, which typically is used as an approximate measure to handle a nonconvexity in an all but convex model, aiming to derive a convex model.


And finally, you can simplify the MATLAB code significantly by vectorizing the code (will not give you a better model though, but improves readability). For instance, summation_of_norm  = sum((y-y_est).^2)/N and similiar for all other stuff

[Mark L. Stone]

norm_nuc is a CVX function.

[Johan Löfberg]

objective = norm_nuc(value(X))+ summation_of_norm ;    

is equivalent to

objective = norm_nuc(someconstantmatrix)+ summation_of_norm ;    

so you are only minimizing summation_of_norm (which is simple as all the nasty nonlinear terms thus mistakingly are omitted)

value simply returns the currently saved value internally in YALMIP (the value you obtain after having done an assignment, or the value after an optimization has been performed.)

[Bhagyashri Telsang]

I have made some changes in the way I was implementing the problem. Now I have removed all the nonlinear terms and expressed it as a linear problem.

However, when I check the diagnostics, it says that the path is Infeasible. Following is the changed code. Kindly go through it and let me know if you can find any flaws. 

Thank you in advance.

n = 3; N = 45; s = 10;

A = randn(n,n);

B = randn(n,1);

C = randn(1,n);

D = randn(1,1);

K = randn(n,1);

u = randn(N,1);

e = randn(N,1);

x = [0; 0; 0];

for k=1:N

    

    y(k) = C*x + D*u(k) + e(k);

    x = (A-K*C)*x + (B-K*D)*u(k) + K*y(k);

    

end

y = y';

%construct Hankel matrix for input and output sequences

U = []; %specify size

for i = 1:s

    for j = 1:N-s+1

        U(i,j) = u(j+i-1);

    end

end

Y = [];

for i = 1:s

    for j = 1:N-s+1

        Y(i,j) = y(j+i-1);

    end

end

v = sdpvar(s,1);

T_u_s = sdpvar(s,s);

for i = 1:s

    for j = 1:s

        if j<=i

            T_u_s(i,j) = v(i-j+1);

        else

            T_u_s(i,j) = 0;

        end

    end

end

w = sdpvar(s-1,1);

T_y_s = sdpvar(s,s);

for i = 1:s

    for j = 1:s

        if j== i

            T_y_s(i,j) = 0;

        elseif j<i

            T_y_s(i,j) = w(i-j);

        else

            T_y_s(i,j) = 0;

        end

    end

end

y_est = sdpvar(N,1);

Y_est = sdpvar(s,(N-s+1));

for i = 1:s

    for j = 1:N-s+1

        Y_est(i,j) = y_est(j+i-1);

    end

end

summation_of_norm = 0;

for i = 1:N

    summation_of_norm = summation_of_norm + (y(i)-y_est(i))^2;

end

summation_of_norm = summation_of_norm/N;

X = Y_est - T_u_s * U - T_y_s * Y;

t1 = norm(X,'nuclear');

t2 = summation_of_norm;

t = t1 + t2;

constraint = [];

objective = t ;    

DIAGNOSTIC = solvesdp(constraint,objective,sdpsettings('solver','sedumi')) ;

    

[Johan Löfberg]

Since your random models typically are unstable, the values in x and y are huge, and then you use the square of y making some values hugehuge, and all solvers I tried choked on the problem due to horrible numerics (infeasible is simply a guess from the solver when the numerics detoriate. Mosek explicitly says it doesn't like the huge values). Hence, you cannot test the model on random numerically unstable examples

Having said that, your code is unnecessarily complex. Some examples of simplifications

T_u_s = flipud((hankel(v(end:-1:1))));

Y_est = hankel(y_est,y_est);
Y_est = temp(1:s,1:N-s+1)

summation_of_norm = (sum(y(1:N)-y_est(1:N)).^2)/N;

and similiar for creating U and Y