Simultaneously inverting many matrices

Dear all, I have many 2-by-2 matrices (which are covariance matrices). I want to invert them all. I'm curious if there's an efficient way of doing this. I thought, maybe, you create a cell, in which each element is one of these matrices, and then use cellfun() in some way to do it. Quintessentially, my question is, is there a way of simultaneously inverting many matrices? I'd appreciate any and all comments. Thank you very much in advance!

Best, John

Answer 1

[Walter Roberson]

2 x 2 you might as well use the formula

D = A(1, 2, :) .* A(2, 1, :) - A(1, 1, :) .* A(2, 2, :);
V11 = -A(2, 2, :) ./ D;
V12 = A(1, 2, :) ./ D;
V21 = A(2, 1, :) ./ D;
V22 = -A(1, 1, :) ./ D;

invs = [V11, V12; V21, V22];

Answer 2

[Tohru Kikawada]

You can leverage Symbolic Math Toolbox to vectorize the calculation.

% Define size of matrices
M=2;
N=10000;
A=rand(M,M,N);

% Calculate the inverse matrices in a loop
invA_loop = zeros(size(A));
tic
for k = 1:N
    invA_loop(:,:,k) = inv(A(:,:,k));
end
disp('Elapsed time in calculation in a loop:');
toc

% Calculate the inverse matrices in a vectorization
As = sym('a', [M,M]); % Define an MxM matrix as a symbolic variable
invAs = reshape(inv(As),[],1); % Solve inverse matrix in symbol
invAfh = matlabFunction(invAs,'Vars',As); % Convert the symbolic function to an anonymous function.
tic
invA_sym = reshape(invAfh(A(1,1,:),A(2,1,:),A(1,2,:),A(2,2,:)),M,M,N);
disp('Elapsed time in the vectorized calculation:');
toc

% Max difference between the results in the loop and the vectorization
disp('Max difference between the elements of the results:');
disp(max(abs(invA_loop(:)-invA_sym(:))))

Results:

Elapsed time in calculation in a loop:

Elapsed time is 0.093938 seconds

Elapsed time in the vectorized calculation:

Elapsed time is 0.001396 seconds

Max difference between the elements of the results:

3.0323e-09

[Lucas Azevedo]

Could you please tell me why there is a different between the elements of the results?

[Tohru Kikawada]

Because function inv performs an LU decomposition which is different from the vectorized one. Please see the documentation for more details.

[Lucas Azevedo]

Humm, I see. By the way, I am using your vectorization and now my code is 2x faster! Thank you!

[Alec Jacobson]

%invA_sym = reshape(invAfh(A(1,1,:),A(2,1,:),A(1,2,:),A(2,2,:)),M,M,N);
Acell = reshape(num2cell(A,3),1,[]);
invA_sym = reshape(invAfh(Acell{:}),M,M,N);

so the code above works for M≠2

[Tohru Kikawada]

Alec, this is great! Thanks for your extension!

Answer 3

[James Tursa]

You might look at this FEX submission by Bruno Luong for solving 2x2 or 3x3 systems:

http://www.mathworks.com/matlabcentral/fileexchange/27762-small-size-linear-solver

Answer 4

[Hugo]

In my experience, using cells is rather slow. Since your matrices are 2x2, then you could simple arrange them in a 3D array, with the first dimension representing the index of each matrix. Let's call this matrix M, which will be of size Nx2x2, N denoting the number of matrices you want to invert.

Now recall that the inverse of a matrix A=[A11,A12;A21,A22] can be computed as

[A22, -A12; -A21, A11] /(A11*A22-A12*A21)

You can implement this for all matrices as follows:

Minv = reshape([M(:,4),-M(:,2),-M(:,3),M(:,1)]./repmat(M(:,1).*M(:,4)-M(:,2).*M(:,3),1,4),N,2,2);

Hope this helps

Hugo

Answer 5

[Azzi Abdelmalek]

Using cellfun will not do it simultaneously. The for loop can be faster. But if you have a Parallel Computing Toolbox, you can do it with parfor