How to approach a optimization problem from experimental data

Hello everyone,

I have acquired values for matrix C(~5000, ~20), through experiments. I believe that the data that I have acquired can be explained by matrix multiplication. In other words A*B=C. Size of A(~5000, ~5) and B(~5,~20). I would like to find the values for the elements in A and B that “best” could explain C. That is the main goal.

I formulated this as an optimization problem: Min sum(sum(((A*B) -C).^2)). There is some linear constraints on elements in B and lower and upper bounds for elements in A and B. I have used the sym toolbox to generate jacobian and hessian for a smaller problem and it works fine. But for a problem of this size the computational time is long. Even when calculating jacobian and hessian using multiple cores. And calculating jacobian and hessian in this way may not be the best option as the hessian is kind of sparse? I am not sure what algorithm would be the best to solve this problem (have used interior-point this far). So suggestion about that is more than welcome. I have used GlobalSearch to try to find the “best” solution.

Is there a better way of finding the values for A and B that explains C? If there is any better way of doing this (objective function, algorithm, you name it) that is less computationally heavy or gives better results, please tell me.

Answer 1

[Matt J]

Maybe you should show what the sym toolbox gave you for the real problem size. But the expressions for the gradient and Hessian should have a really simple matrix form for this problem - easy to develop by hand without the sym toolbox and very quick to evaluate. For example, the gradients with respect to A and B are

    gradient_wrt_A = (A*B-C)*B.';
    gradient_wrt_B = A.'(A*B-C);

[Niklas handin]

Dear Matt,

Thanks for your replay. I have now added a short example of my code. I did decrees the size of the matrixes so it wouldn’t take the whole day to run it. Thanks for the great suggestion about the gradients. But isn’t it missing a factor 2?

gradient_wrt_A = 2*((A*B-C)*B.');
gradient_wrt_B = 2*(A.')*(A*B-C);

I guess it’s time for me to refresh gradient and hessian calculations. If you have any other suggestions to improve/change my approach it would very appreciated.

Best Niklas

[Matt J]

But isn’t it missing a factor 2?

Yes, I guess so, although it would be more traditional for you to write your objective function as

   f(A,B)=0.5*norm(A*B-C,'fro')^2,

in which case the factor of 2 in the gradient goes away.