Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Multiple backslash on the same rhs

Subject: Multiple backslash on the same rhs

From: Bruno Luong

Date: 14 Sep, 2008 13:49:01

Message: 1 of 7

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

Subject: Multiple backslash on the same rhs

From: John D'Errico

Date: 14 Sep, 2008 14:50:03

Message: 2 of 7

"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

Subject: Multiple backslash on the same rhs

From: Tim Davis

Date: 14 Sep, 2008 19:11:01

Message: 3 of 7

"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).

Subject: Multiple backslash on the same rhs

From: Bruno Luong

Date: 14 Sep, 2008 22:05:04

Message: 4 of 7

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

Subject: Multiple backslash on the same rhs

From: Bruno Luong

Date: 14 Sep, 2008 22:33:02

Message: 5 of 7


>
>
> 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

Subject: Multiple backslash on the same rhs

From: Tim Davis

Date: 15 Sep, 2008 00:57:02

Message: 6 of 7

"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.

Subject: Multiple backslash on the same rhs

From: Bruno Luong

Date: 15 Sep, 2008 05:40:04

Message: 7 of 7

"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

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us