Thread Subject:

Multiple backslash on the same rhs

I have p (m x n)-matrices ranging in a 3D arrays M(1:m,1:n,1:p)

I would like to apply multiple backslash operators of these matrices on the same rhs, then ranging the result in a 2D array X(1:n,1:p).

In my application typically m and n are small (<=5) but p can be quite large (hundreds to few thousands).

So far I use loop:

m = 3;
n = 2;
p = 5;

% Test matrix and Data
M = rand(m,n,p);
rhs = rand(m,1);

% Solving
X = zeros(n,p);
for k=1:P
  X(:,k) = M(:,:,k) \ rhs;
end

Is there any better (faster) way?

This kind of problem seems to be a recurrent one for me, but I couldn't find any better way to accomplish it.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gaj4od$t5i$1@fred.mathworks.com>...
> I have p (m x n)-matrices ranging in a 3D arrays M(1:m,1:n,1:p)
>
> I would like to apply multiple backslash operators of these matrices on the same rhs, then ranging the result in a 2D array X(1:n,1:p).
>
> In my application typically m and n are small (<=5) but p can be quite large (hundreds to few thousands).
>
> So far I use loop:
>
> m = 3;
> n = 2;
> p = 5;
>
> % Test matrix and Data
> M = rand(m,n,p);
> rhs = rand(m,1);
>
> % Solving
> X = zeros(n,p);
> for k=1:P
> X(:,k) = M(:,:,k) \ rhs;
> end
>
> Is there any better (faster) way?
>
> This kind of problem seems to be a recurrent one for me, but I couldn't find any better way to accomplish it.
>
> Bruno

Probably not. Certainly not unless the matrices
are strongly related to each other. For example,
if M(:,:,i) and M(:,:,i+1) are the same, but for a
rank 1 update of some form.

Look into qrupdate and qrdelete. If they will not
help you, then I see no solution.

John

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gaj4od$t5i$1@fred.mathworks.com>...
> I have p (m x n)-matrices ranging in a 3D arrays M(1:m,1:n,1:p)
>
> I would like to apply multiple backslash operators of these matrices on the same rhs, then ranging the result in a 2D array X(1:n,1:p).
>
> In my application typically m and n are small (<=5) but p can be quite large (hundreds to few thousands).
>
> So far I use loop:
>
> m = 3;
> n = 2;
> p = 5;
>
> % Test matrix and Data
> M = rand(m,n,p);
> rhs = rand(m,1);
>
> % Solving
> X = zeros(n,p);
> for k=1:P
> X(:,k) = M(:,:,k) \ rhs;
> end
>
> Is there any better (faster) way?
>
> This kind of problem seems to be a recurrent one for me, but I couldn't find any better way to accomplish it.
>
> Bruno

Try creating a block diagonal matrix, with p diagonal blocks each of size m-by-n. To do that, first create a list of nonzeros in 3 vectors of length s-by-1 where s=m*n*p (row indices, column indices, and values). The last one is just M(:). Then do

A = sparse (I,J,M(:),m*p,n*p);

Then do

X = A\rhs(:);
X = reshape (X, something here to convert a vector to a matrix).

Thank you Tim and John,

Here is a test with Tim's suggestion

m = 5;
n = 3;
p = 1e4;

% For loop
tic
M = rand(m,n,p);
rhs = rand(m,1);
X1=zeros(n,p);
for k=1:p
  X1(:,k) = M(:,:,k) \ rhs;
end
toc % 0.228802 seconds.

% Sparse matrix
%[m n p]=size(M);
tic
I = repmat(reshape(1:m*p,m,1,p),[1 n 1]);
J = repmat(reshape(1:n*p,1,n,p),[m 1 1]);
A = sparse(I(:),J(:),M(:));
X2 = reshape(A \ repmat(rhs,p,1), n, p);
toc % 0.158009 seconds.


It seems sparse is faster!!!

Now I have a crazy idea: can't we just implement the inverse with determinant formula and apply the same formula for all matrices?

Is there any Matlab tool for implementing inverse using determinant (I know it must be very inefficient method for large matrix, but I would think it works OK on matrix with dimension less than 5).

Bruno

>
>
> It seems sparse is faster!!!
>

Sorry the first "tic" should move down, right bebore X1=...

The conclusion is still the same: sparse is faster than for-loop

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gak1qg$4d8$1@fred.mathworks.com>...
> Thank you Tim and John,
>
> Here is a test with Tim's suggestion
>
> m = 5;
> n = 3;
> p = 1e4;
>
> % For loop
> tic
> M = rand(m,n,p);
> rhs = rand(m,1);
> X1=zeros(n,p);
> for k=1:p
> X1(:,k) = M(:,:,k) \ rhs;
> end
> toc % 0.228802 seconds.
>
> % Sparse matrix
> %[m n p]=size(M);
> tic
> I = repmat(reshape(1:m*p,m,1,p),[1 n 1]);
> J = repmat(reshape(1:n*p,1,n,p),[m 1 1]);
> A = sparse(I(:),J(:),M(:));
> X2 = reshape(A \ repmat(rhs,p,1), n, p);
> toc % 0.158009 seconds.
>
>
> It seems sparse is faster!!!

Maybe not. You charged the time to compute M and rhs to the for loop, but not to the sparse version.

> Now I have a crazy idea: can't we just implement the inverse with determinant formula and apply the same formula for all matrices?

I would not recommend it. Besides, the matrix is rectangular; I don't think the determinant is defined for rectangular matrices.

If you want to repeat this with the same A but different right-hand-sides, do:

[Q,R]=qr(A) ;

once, and then do

x = R\(Q'*b)

repeatedly, where A is your sparse matrix above.

Normally, Q has way too many nonzeros for it to be efficient, but for this block diagonal matrix A, Q is pretty sparse.

"Tim Davis" <davis@cise.ufl.edu> wrote in message <gakbsu$c7k$1@fred.mathworks.com>...

> I would not recommend it. Besides, the matrix is rectangular; I don't think the determinant is defined for rectangular matrices.

What I have in mind is apply inverse on M(:,:,k).'*M(:,:,k), As you know, this comes from the normal equation. My system is underdeermined one, so it is the same as "\".

Bruno