Fast maximum likelihood matrix decomposition

I need to decompose a matrix (H) into two matrices (a0 and a0'), so that H=a0*a0'. There is no analytical solution to the problem, so I'm using a maximum likelihood estimation. More precisely, I'm using the following function:

    function [a0best]=decomp(H);

    %constraints for the matrix a0:
    lb=[-inf 0 0 ; -inf -inf 0; -inf -inf -inf];
    ub=[inf 0 0 ; inf inf 0; inf inf inf]; 

    %options for the solver:
    options=optimset('Algorithm','sqp','Display','off');

    %starting value for the matrix a0:
    a0=tril(rand(3));

    %minimization:
    a0best=fmincon(@llfazero,a0,[],[],[],[],lb,ub,[],options);
    a0best=eye(3)/a0best;    

    %the likelihood function:
    function [llf]=llfazero(a0);   
    llf = -log(abs(det(a0))) + 0.5* sum(diag(a0*H*a0'));   
    end
    end

The function works fine, but very slowly. According to my estimations, it would take me 5 days to run the whole code with this function. So, is there a faster way to run this minimization? I stress that there is no analytical solution to my problem, unlike in the above (simplified) code, where the solution is, obviously, the Cholesky decomposition. Many thanks.

Answer 1

[Matt J]

According to my estimations, it would take me 5 days to run the whole code with this function.

You should elaborate on how you estimate that and what "slow" means in this case. I assume you're seeing slow convergence and not long computation time per iteration. If the space of unknowns is really 3x3 matrices a0, that objective function you've shown should be trivial to compute quickly.

A few remarks, though

1. It is a bad idea to try to force some of your a0(i,j) to be zero using lb and ub bounds. If some of the a0(i,j) have known values, you should just exclude them from the list of unknowns, e.g., like in the revised version of the code below.
2. It would probably be better to add a positivity constraint on det(a0) and get rid of abs(det(a0)) in your objective function. Otherwise, the formulation of your problem is not differentiable.

      pbest=fmincon(@llfazero,rand(6,1),[],[],[],[],[],[],nonlcon,options);

      %the likelihood function:
      function [llf]=llfazero(p);   

        a0(tril(true(3)))=p; 
            a0=reshape(a0,3,3);

        llf = -log(p(1)*p(4)*p(6)) + 0.5* sum(diag(a0*H*a0'));   
      end

      function [c,ceq]=nonlcon(p)

       ceq=[];
       c=-p(1)*p(4)*p(6); %positive determinant

      end

[bran]

Thanks, Matt, your suggestions reduced my computing time by nearly 20%. The function is simple, sure, but I am running some simulation and I need to compute it million times (or maybe more). The rest of the code should be fine - when I run it with analytical decomposition of the matrix, a number of iterations takes 4 seconds, but when I run it with maximum likelihood decomposition, it takes 40 seconds. So, the slow part is this decomposition, not the rest of the code.