Matlab taking too much time to give output

I'm solving a convex optimization problem which has to minimize a variable matrix of 498x498 size. As I run it, it is taking forever to give final value of x. When I terminate the process in between it shows like this

my code is

function f= coheren(x)

   load ('Dictionary.mat','D');
   Z = (D)*(D');
   [V, U,W]=svd(Z);
   a = mtimes(U,V');
   b = mtimes(V,U);
   f = (norm(U-(a*(x')*(x)*b)))^2;
  end

function [c,ceq] = cohcon(x)

   load ('Dictionary.mat','D');
   [m,n]= size(D);%size(x)
   k = D*(D');
   T = (x)*k*(x');
   [V1,U1,W1] =svd(T);
   ceq = rank(U1)-m;
   c= [];
  end

Answer 1

[Walter Roberson]

Never load() inside a cost function. Load the data outside the function and pass the data in.

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

In your case, for coheren(), you could pre-calculate U, a, and b, and pass those all in. For cohcon, you could pre-calculate k = D*(D') and pass that in.

[Pushkar Khatri]

But, it still not helping. I followed as you suggested but it is still taking too long

functions are -

function f= coheren(x,a,b,U)

   f = (norm(U-(a*(x')*(x)*b)))^2;

end

function [c,ceq] = cohcon(x,m,Z)

   T = (x)*Z*(x');
   [~,U1,~] =svd(T);
   c = rank(U1)-m;
   ceq= [];

end

[Walter Roberson]

Can you attach representative data for us to test with?

... And remember, we cannot paste images of code into MATLAB for execution.

[Pushkar Khatri]

Here is the dictionary library having a D (498x221) matrix. Command line codes:

 >>load('Dictionary.mat','D');
 Z = (D)*(D');
 [V, U,W]=svd(Z);
 a = mtimes(U,V');
 b = mtimes(V,U);
 clear W;
 [m,~] = size(D);
>>  options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
>> fun = @(x)coheren(x,a,b,U);
>> nonlcon = @(x)cohcon(x,m,Z);
>> x = fmincon(fun,zeros(498),[],[],[],[],[],[],nonlcon,options)

[Walter Roberson]

The initial gradient and Hessian estimation takes a long time!

[Pushkar Khatri]

So what should I do now and how long should stay awake for the output? Will it ever give the result?

[Walter Roberson]

>> cohdrive
 Iter  Func-count            Fval   Feasibility   Step Length       Norm of   First-order  
                                                                       step    optimality
    0      248005    1.244865e+09     0.000e+00     1.000e+00     0.000e+00     7.856e+03  
Error using sqpLineSearchMex
Requested 248004x248004 (458.3GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more
information.
Error in sqpInterface
Error in fmincon (line 823)
    [X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = sqpInterface(funfcn,X,full(A),full(B),full(Aeq),full(Beq), ...
Error in cohdrive (line 11)
 x = fmincon(fun,zeros(498),[],[],[],[],[],[],nonlcon,options)

[Walter Roberson]

The same thing happened when I switched to interior-point, except that it did not get to displaying even one iteration.

[Pushkar Khatri]

Yes,for normal or smaller datasets ,the sqp algorithm is faster than interior-point algorithm. And gives the final x in lesser no. of iterations.

[Walter Roberson]

Do you have a comment on the analysis I did for the nonlinear constraints?

[Pushkar Khatri]

First of all, thanks a lot for your help!! The rank of the dictionary,used here, is 221 and after svd the diagonal matrix(here matrix U) formed is not a full rank matrix. Guessing in the analogous way,I also think <=0 constraint will allways be true for x. But fmincon requires at least one constraint to carry on the convex optimization. Without that case, hour can it possibly operate

[Walter Roberson]

Convex optimization would be more associated with linear constraints rather than with nonlinear constraints.

Answer 2

[Aakash Deep]

Hey Pushkar,

Solving a convex optimization problem is a big task. Yes, it can take a long time to converge. You can not minimize the time of the algorithm but what you can do is, you can change your converging threshold.

As you know, the convergence plot of convex optimization problem decrease exponentially and hence the major components are already gets covered in the starting. So, you can give a threshold (try to keep the threshold at or below 0.1) such that whenever your error value gets below the threshold then stop the algorithm.

[Pushkar Khatri]

I didn't understand completely. Can you please explain more about converging threshold. Are you saying about the starting point that I took as zeros(498)?

[Aakash Deep]

Okay, no problem. There is another simple approach for this. Try to use maximum iteration parameter in your optimoptions to stop your algorithm (50 iterations are more than enough to recover a signal efficiently but obviously you can change it according to your need).

For further details, refer this,

https://www.mathworks.com/help/optim/ug/optimoptions.html#btm5a5a-4

Hope this helps :)

[Walter Roberson]

It is not possible to complete even one iteration, because of the (498*498) x (498*498) matrix that is needed for hessian estimation. The interaction of every variable with every other variable must be estimated.

I wonder if this could be rewritten as "large-scale" ?

[Aakash Deep]

I don't think the size of the matrix will create any problem (50 iteration was just a suggested value, obviously it will be different for everyone). By this, the major problem of code not stopping will get solve.

[Walter Roberson]

Aakash Deep, I actually executed the code. Before the second iteration ('sqp') or before the first iteration ('interior-point') the error message is given,

Requested 248004x248004 (458.3GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more
information.

248004 is exactly 498^2, which is the number of variables involved because zeros(498) is passed as the x value to fmincon. Hessian estimation requires estimating the effect of each variable with respect to each other variable, which would lead to a 248004 x 248004 matrix, which is 458.3GB . I know that I do not have that much memory!

[Aakash Deep]

Yes, I agree with your point. So in that case, is there any solution for this problem?

[Walter Roberson]

The only hope I see while using fmincon() is if the algorithm can be expressed for use in large scale solvers.

When I look at the nonlinear constraint, after some thinking, it looks to me as if it testing whether the matrix is singular in the sense of having at least one zero singular value. It calculates the rank of the 498 x 498 diagonal matrix, with the rank being the same as the number of non-zero diagonal elements, and it subtracts off the number of rows of the D matrix (498) and returns that. That would be negative in the case where there was at least one zero singular value, and would be 0 in the cast where the passed in matrix had no non-zero singular values.

... But -- nonlinear inequalities are consider satisfied if the output is <= 0, which is always going to be the case because rank minus 498 with rank maximum 498, can never be positive. With the two cases distinguished in theory seeming to be either singular or not, then I suspect that 1 should have been added to the result of the test so that the result would be positive in the case where the input was full rank.

When I look at https://www.mathworks.com/help/matlab/ref/svd.html#bu2_4_9 it says you can calculate the rank of a matrix by calculating the number of non-zero singular values -- which is what appears to be happening here. But it also says that rank() is provided for that purpose. This suggest that even if the nonlinear constraint is still required, that the SVD calculation followed by finding the rank of the diagonal, could be replaced by simply calling rank() on the matrix.

As I look through the way that matrix is being constructed, rank(Z) is 215, and I do not think it is possible for the calculations it is doing to create a greater rank matrix. So my thought is that at least for this D matrix, that the nonlinear constraint will always return negative (meaning constraint not violated). If my logic is correct then the nonlinear constraint call could be removed. That could be important for the purpose of converting to large scale algorithms.

Unfortunately at the moment I have no idea how one would calculate the hessian pattern for the large scale algorithms.

Answer 3

[Rishabh Rathore]

For the convergence problem, the upper bound could be used on the number of iterations, as suggested by Aakash. As for the problem with size of array being too large, Tall arrays could be used to overcome the issue as Tall arrays are not loaded in the memory at once.

You could refer:- Tall Arrays

Note:- the computation would take time since the size of the data to be manipulated and the number of calculations to be performed is very high.