How can I optimise solving a linear system for multiple points.

I want to solve a linear problem:

N=30;
%vector x is filled from a file and has size 2,30
%vector y is filled from a file and is a multidimensional array with size 100,100,3,30
fixedpoint1 = 50;
fixedpoint2 = 45;
for i=1:N
   Matrix(i,:) = [x(i,1)^2, x(i,2)^2, x(i,1), x(i,2), 1]
   z(i) = y(fixedpoint1, fixedpoint2, 1, i)
end

solution=Matrix\z'

I would like to do this for every fixedpoint1 and fixedpoint2 (100x100 times) and I suspect that matlab has some powerful shortcuts for doing so, that are not just creating two for loops. Do you have any suggestions?

I know that I can solve this problem also using SVD, would it be better? I cannot see how it could be.

Any suggestion to make the code more "matlab-like" is welcome.

Answer 1

[Cam Salzberger]

Hello Andrea,

I understand that you are computing some least-squares solutions for a multi-dimensional array, and would like to know how to do it quickly. I have two examples, where I compute the "solution" values for all of the "fixedpoints". The "solution" vector for each combination of the two fixed points is stored in a three-dimensional array, such that:

solutionArr(fixedpoint1,fixedpoint2,:) = the "solution" from your code

Firstly, using vectorized code to create the inputs to the least-squares solver (the "\" operator ) will remove the "for" loop that you are currently showing in your code. The use of the "permute" function is quite useful in this regard. This can be seen in the code below:

Matrix = [x.^2 x ones(size(x,1),1)];
solutionArr = zeros(size(y,1),size(y,2),size(Matrix,2));
for r = 1:size(y,1)
    for c = 1:size(y,2)
        z = permute(y(r,c,1,:),[4 2 1 3]);
        solutionArr(r,c,:) = permute(Matrix\z,[3 2 1]);
    end
end

To improve the execution time of the code even more, you can take advantage of the ability to use the "x = A\b" least-squares solution, but when "x" and "b" are matrices. Removing another "for" loop will speed up the code even more.

Matrix = [x.^2 x ones(size(x,1),1)];
solutionArr = zeros(size(y,1),size(y,2),size(Matrix,2));
for r = 1:size(y,1)
    Z = permute(y(r,:,1,:),[4 2 1 3]);
    solutionArr(r,:,:) = permute(Matrix\Z,[3 2 1]);
end

You may also wish to look into orthogonal-triangular decomposition or LU matrix factorization . Both of these may calculate the least-square solution in less time than the "\" operator.

I hope that this helps.

-Cam