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:
overdetermenied system

Subject: overdetermenied system

From: Awusi Kavuma

Date: 13 Jan, 2009 11:45:04

Message: 1 of 12

Hello there

Can someone help with this overdetermenied system

M1(x,y) = A(x,y).*M11(x,y) + B(x,y).*M10(x,y)
M2(x,y) = A(x,y).*M22(x,y) + B(x,y).*M20(x,y)
M3(x,y) = A(x,y).*M33(x,y) + B(x,y).*M30(x,y)
M4(x,y) = A(x,y).*M44(x,y) + B(x,y).*M40(x,y)

All Mi(x,y) are 512x384 matrices, real and known.

The problem is linear in A and B and for specific points the least square using backslash technique can be used to solve for the coefficients. My problem is to extend the method to find coefficients A(x,y) and B(x,y) for entire system.

Cheers!

Kavuma

Subject: overdetermenied system

From: John D'Errico

Date: 13 Jan, 2009 12:53:02

Message: 2 of 12

"Awusi Kavuma" <kavumawusi@yahoo.com> wrote in message <gkhurv$nv1$1@fred.mathworks.com>...
> Hello there
>
> Can someone help with this overdetermenied system
>
> M1(x,y) = A(x,y).*M11(x,y) + B(x,y).*M10(x,y)
> M2(x,y) = A(x,y).*M22(x,y) + B(x,y).*M20(x,y)
> M3(x,y) = A(x,y).*M33(x,y) + B(x,y).*M30(x,y)
> M4(x,y) = A(x,y).*M44(x,y) + B(x,y).*M40(x,y)
>
> All Mi(x,y) are 512x384 matrices, real and known.
>
> The problem is linear in A and B and for specific points the least square using backslash technique can be used to solve for the coefficients. My problem is to extend the method to find coefficients A(x,y) and B(x,y) for entire system.
>

This reduces to a linear block diagonal system
with 4x2 blocks. Backslash should be able to
solve it efficiently.

John

Subject: overdetermenied system

From: Roger Stafford

Date: 13 Jan, 2009 13:22:02

Message: 3 of 12

"Awusi Kavuma" <kavumawusi@yahoo.com> wrote in message <gkhurv$nv1$1@fred.mathworks.com>...
> Hello there
>
> Can someone help with this overdetermenied system
>
> M1(x,y) = A(x,y).*M11(x,y) + B(x,y).*M10(x,y)
> M2(x,y) = A(x,y).*M22(x,y) + B(x,y).*M20(x,y)
> M3(x,y) = A(x,y).*M33(x,y) + B(x,y).*M30(x,y)
> M4(x,y) = A(x,y).*M44(x,y) + B(x,y).*M40(x,y)
>
> All Mi(x,y) are 512x384 matrices, real and known.
>
> The problem is linear in A and B and for specific points the least square using backslash technique can be used to solve for the coefficients. My problem is to extend the method to find coefficients A(x,y) and B(x,y) for entire system.
>
> Cheers!
>
> Kavuma

  Kavuma, the way you have written your equations using element-wise multiplication, ".*", would imply that A and B are also each 512 x 384 matrices. Is that what you intend? When you say "for specific points", do you mean for each of the specific elements among the total of 512*384 = 196608 elements? Is that the meaning the notation Mi(x,y), A(x,y), etc. is meant to convey?

  If that is the case, what you have is a system of 196608 instances of four equations and two unknowns, each of which you want solved by least squares. For each such individual instance you could use the backslash operator. However I doubt if this operator can be made to work for all elements simultaneously in a single expression. Presumably either you would have to do it within nested for-loops for each individual element, or else write out the specific solution to the four-linear-equations-with-two-unknowns problem as a single vectorized expression.

  With so few unknowns and equations the latter isn't impossibly difficult to do, but before going on with this, it would be advisable for you to tell us whether the above interpretation is correct.

Roger Stafford

Subject: overdetermenied system

From: John D'Errico

Date: 13 Jan, 2009 13:53:02

Message: 4 of 12

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gki4hq$1b9$1@fred.mathworks.com>...

> Kavuma, the way you have written your equations using element-wise multiplication, ".*", would imply that A and B are also each 512 x 384 matrices. Is that what you intend? When you say "for specific points", do you mean for each of the specific elements among the total of 512*384 = 196608 elements? Is that the meaning the notation Mi(x,y), A(x,y), etc. is meant to convey?
>
> If that is the case, what you have is a system of 196608 instances of four equations and two unknowns, each of which you want solved by least squares. For each such individual instance you could use the backslash operator. However I doubt if this operator can be made to work for all elements simultaneously in a single expression.

Actually, if you build it as a sparse block diagonal
matrix, I think that backslash might survive the
ordeal. At the very least, you could do it in chunks,
with perhaps 512 blocks in each chunk. The matrix
problem MUST be done as a sparse one however.

John

Subject: overdetermenied system

From: Awusi Kavuma

Date: 13 Jan, 2009 14:22:01

Message: 5 of 12

Hi Rogers

1) Yes each of A and B should also be 512x384.
2) True each of the A(x,y) and M(x,y) has 512x384 = 196608
3) For specific point I mean forexample the cental point is (256,192). I can extract the reading at each of the matrices at 256,192 and form 4 linear equations and use the back slash to find the coefficients. My problem arise when it comes to generalise it to find all the coefficients.
4) There are 4 equations with two unknown matrices.

Let me hope this clarifies your questions and thank again for more of your coments.


kavuma

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gki4hq$1b9$1@fred.mathworks.com>...
> "Awusi Kavuma" <kavumawusi@yahoo.com> wrote in message <gkhurv$nv1$1@fred.mathworks.com>...
> > Hello there
> >
> > Can someone help with this overdetermenied system
> >
> > M1(x,y) = A(x,y).*M11(x,y) + B(x,y).*M10(x,y)
> > M2(x,y) = A(x,y).*M22(x,y) + B(x,y).*M20(x,y)
> > M3(x,y) = A(x,y).*M33(x,y) + B(x,y).*M30(x,y)
> > M4(x,y) = A(x,y).*M44(x,y) + B(x,y).*M40(x,y)
> >
> > All Mi(x,y) are 512x384 matrices, real and known.
> >
> > The problem is linear in A and B and for specific points the least square using backslash technique can be used to solve for the coefficients. My problem is to extend the method to find coefficients A(x,y) and B(x,y) for entire system.
> >
> > Cheers!
> >
> > Kavuma
>
> Kavuma, the way you have written your equations using element-wise multiplication, ".*", would imply that A and B are also each 512 x 384 matrices. Is that what you intend? When you say "for specific points", do you mean for each of the specific elements among the total of 512*384 = 196608 elements? Is that the meaning the notation Mi(x,y), A(x,y), etc. is meant to convey?
>
> If that is the case, what you have is a system of 196608 instances of four equations and two unknowns, each of which you want solved by least squares. For each such individual instance you could use the backslash operator. However I doubt if this operator can be made to work for all elements simultaneously in a single expression. Presumably either you would have to do it within nested for-loops for each individual element, or else write out the specific solution to the four-linear-equations-with-two-unknowns problem as a single vectorized expression.
>
> With so few unknowns and equations the latter isn't impossibly difficult to do, but before going on with this, it would be advisable for you to tell us whether the above interpretation is correct.
>
> Roger Stafford

Subject: overdetermenied system

From: Roger Stafford

Date: 13 Jan, 2009 14:26:02

Message: 6 of 12

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gki6bu$6j$1@fred.mathworks.com>...
> Actually, if you build it as a sparse block diagonal
> matrix, I think that backslash might survive the
> ordeal. At the very least, you could do it in chunks,
> with perhaps 512 blocks in each chunk. The matrix
> problem MUST be done as a sparse one however.
>
> John

  Yes, I expect you are right, John. It could be done that way. I shouldn't make such generalized assertions.

  I'm not so sure about its efficiency, though, as opposed to 512*384 individual 4x2 backslashes. The single block diagonal backslash would have to be smart enough not to waste significant overhead time on all the implied zero "off-diagonals". The solution to a single 4x2 problem should be very fast, especially if one writes it out as an explicit vectorized expression in terms of the four basic arithmetic operations without using matlab's backslash operator itself.

Roger Stafford

Subject: overdetermenied system

From: Bruno Luong

Date: 13 Jan, 2009 14:40:20

Message: 7 of 12

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gki89q$85e$1@fred.mathworks.com>...

>
> I'm not so sure about its efficiency, though, as opposed to 512*384 individual 4x2 backslashes. The single block diagonal backslash would have to be smart enough not to waste significant overhead time on all the implied zero "off-diagonals". The solution to a single 4x2 problem should be very fast, especially if one writes it out as an explicit vectorized expression in terms of the four basic arithmetic operations without using matlab's backslash operator itself.

Roger,

Yes, "\" is smart enough. It is slighly faster than for-loop, see this thread that use the same trick: http://www.mathworks.de/matlabcentral/newsreader/view_thread/235989

Bruno

Subject: overdetermenied system

From: John D'Errico

Date: 13 Jan, 2009 15:03:01

Message: 8 of 12

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gki89q$85e$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gki6bu$6j$1@fred.mathworks.com>...
> > Actually, if you build it as a sparse block diagonal
> > matrix, I think that backslash might survive the
> > ordeal. At the very least, you could do it in chunks,
> > with perhaps 512 blocks in each chunk. The matrix
> > problem MUST be done as a sparse one however.
> >
> > John
>
> Yes, I expect you are right, John. It could be done that way. I shouldn't make such generalized assertions.
>
> I'm not so sure about its efficiency, though, as opposed to 512*384 individual 4x2 backslashes. The single block diagonal backslash would have to be smart enough not to waste significant overhead time on all the implied zero "off-diagonals". The solution to a single 4x2 problem should be very fast, especially if one writes it out as an explicit vectorized expression in terms of the four basic arithmetic operations without using matlab's backslash operator itself.
>
> Roger Stafford

If your matrix is sparse, then \ really is pretty smart.

blocks = rand(4,2,1000);
A = blktridiag(blocks,zeros(4,2,999),zeros(4,2,999));
rhs = rand(4000,1);

tic
coef = zeros(2,1000);
for i = 1:1000
coef(:,i) = blocks(:,:,i)\rhs((1:4) + 4*(i-1));
end
toc
Elapsed time is 0.105407 seconds.

tic,coef = A\rhs;toc
Elapsed time is 0.044134 seconds.

I will admit that the simple loop is pretty fast though,
especially if I add in the time to create the sparse
block diagonal matrix. blktridiag is not as efficient
as it should be here, since it fills in explicit zeros
into the upper and lower blocks.

John

Subject: overdetermenied system

From: Roger Stafford

Date: 13 Jan, 2009 16:12:02

Message: 9 of 12

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gkiaf5$300$1@fred.mathworks.com>...
> If your matrix is sparse, then \ really is pretty smart.
>
> blocks = rand(4,2,1000);
> A = blktridiag(blocks,zeros(4,2,999),zeros(4,2,999));
> rhs = rand(4000,1);
>
> tic
> coef = zeros(2,1000);
> for i = 1:1000
> coef(:,i) = blocks(:,:,i)\rhs((1:4) + 4*(i-1));
> end
> toc
> Elapsed time is 0.105407 seconds.
>
> tic,coef = A\rhs;toc
> Elapsed time is 0.044134 seconds.
>
> I will admit that the simple loop is pretty fast though,
> especially if I add in the time to create the sparse
> block diagonal matrix. blktridiag is not as efficient
> as it should be here, since it fills in explicit zeros
> into the upper and lower blocks.
>
> John

  OK, John and Bruno, I surrender. Your sparse block diagonal method seems to win the day. I was unaware of the 9/14/08 thread mentioned by Bruno.

  Nevertheless, it somehow offends my sense of how such an algorithm ought to work. Instead of a procedure that has to be smart enough to ignore off-diagonal zeros in a sparse matrix, I dream of one that never has to deal with them in the first place. Ideally there ought to be a method of carrying out the operations involved in each independent block in a repetitive manner that avoids the inefficiencies of for-loop indexing and its overhead. I suspect the apparent efficiency here of the single backslash operation as opposed to multiple smaller ones lies in being able to do the rough equivalent of for-looping down at a lower language level, which is therefore more efficient than having to use a matlab compiler-generated for-loop structure. I can see no other reason the block diagonal method could win out over simply executing the blocks singly.

Roger Stafford (The dreamer)

Subject: overdetermenied system

From: John D'Errico

Date: 13 Jan, 2009 16:24:01

Message: 10 of 12

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gkiegi$alv$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gkiaf5$300$1@fred.mathworks.com>...
> > If your matrix is sparse, then \ really is pretty smart.
> >
> > blocks = rand(4,2,1000);
> > A = blktridiag(blocks,zeros(4,2,999),zeros(4,2,999));
> > rhs = rand(4000,1);
> >
> > tic
> > coef = zeros(2,1000);
> > for i = 1:1000
> > coef(:,i) = blocks(:,:,i)\rhs((1:4) + 4*(i-1));
> > end
> > toc
> > Elapsed time is 0.105407 seconds.
> >
> > tic,coef = A\rhs;toc
> > Elapsed time is 0.044134 seconds.
> >
> > I will admit that the simple loop is pretty fast though,
> > especially if I add in the time to create the sparse
> > block diagonal matrix. blktridiag is not as efficient
> > as it should be here, since it fills in explicit zeros
> > into the upper and lower blocks.
> >
> > John
>
> OK, John and Bruno, I surrender. Your sparse block diagonal method seems to win the day. I was unaware of the 9/14/08 thread mentioned by Bruno.
>
> Nevertheless, it somehow offends my sense of how such an algorithm ought to work. Instead of a procedure that has to be smart enough to ignore off-diagonal zeros in a sparse matrix, I dream of one that never has to deal with them in the first place. Ideally there ought to be a method of carrying out the operations involved in each independent block in a repetitive manner that avoids the inefficiencies of for-loop indexing and its overhead. I suspect the apparent efficiency here of the single backslash operation as opposed to multiple smaller ones lies in being able to do the rough equivalent of for-looping down at a lower language level, which is therefore more efficient than having to use a matlab compiler-generated for-loop structure. I can see no other reason the block diagonal method could win out over simply executing the blocks singly.
>
> Roger Stafford (The dreamer)

Yes, you can think of the block diagonal method
as just avoiding explicit loops, trading them for
implicit loops in blas calls.

One could also probably have used cellfun or arrayfun
to solve this problem, but the implicit loops used
there are not at quite as low a level as the block
diagonal solution.

John

Subject: overdetermenied system

From: Awusi Kavuma

Date: 14 Jan, 2009 13:01:04

Message: 11 of 12

Dear John

Problem solved.

Thank you so much for your assistance.
Special thanks to John, Rogers, Bruno and Tim

Kavuma

Given these arrays of data,


M1 = [23 108;105 186], M11 = [5 16; 7 18], M10 = [9 10; 11 12]
M2 = [19 156;90 202], M22 = [17 8;19 2], M20 = [1 22; 2 24]
M3 = [-6 -18;81 100], M33 = [-10 4;8 4], M30 = [2 -5; 7 10]
M4 = [6 -3;55 77], M44 = [-18 7;5 1], M40 = [12 -4; 5 9]


Combine them all into one array.
  Mij = [M11(:),M22(:),M33(:),M44(:),M10(:),M20(:),M30(:),M40(:)];
You have n blocks.
  n = numel(M11);
Reshape into a 3d array.
  Mij = reshape(Mij',[4 2 n]);
Then use my blktridiag to make it into a sparse block diagonal
array that backslash can use.
  A = blktridiag(Mij,zeros(4,2,n-1),zeros(4,2,n-1))
Generate the right hand side for backslash in the same way.
  rhs = [M1(:),M2(:),M3(:),M4(:)]';
  rhs = rhs(:);
And solve using backslash.
  AB = A\rhs
A = reshape(AB(1:2:end),size(M11))
B = reshape(AB(2:2:end),size(M11))



"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gkif71$t56$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gkiegi$alv$1@fred.mathworks.com>...
> > "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gkiaf5$300$1@fred.mathworks.com>...
> > > If your matrix is sparse, then \ really is pretty smart.
> > >
> > > blocks = rand(4,2,1000);
> > > A = blktridiag(blocks,zeros(4,2,999),zeros(4,2,999));
> > > rhs = rand(4000,1);
> > >
> > > tic
> > > coef = zeros(2,1000);
> > > for i = 1:1000
> > > coef(:,i) = blocks(:,:,i)\rhs((1:4) + 4*(i-1));
> > > end
> > > toc
> > > Elapsed time is 0.105407 seconds.
> > >
> > > tic,coef = A\rhs;toc
> > > Elapsed time is 0.044134 seconds.
> > >
> > > I will admit that the simple loop is pretty fast though,
> > > especially if I add in the time to create the sparse
> > > block diagonal matrix. blktridiag is not as efficient
> > > as it should be here, since it fills in explicit zeros
> > > into the upper and lower blocks.
> > >
> > > John
> >
> > OK, John and Bruno, I surrender. Your sparse block diagonal method seems to win the day. I was unaware of the 9/14/08 thread mentioned by Bruno.
> >
> > Nevertheless, it somehow offends my sense of how such an algorithm ought to work. Instead of a procedure that has to be smart enough to ignore off-diagonal zeros in a sparse matrix, I dream of one that never has to deal with them in the first place. Ideally there ought to be a method of carrying out the operations involved in each independent block in a repetitive manner that avoids the inefficiencies of for-loop indexing and its overhead. I suspect the apparent efficiency here of the single backslash operation as opposed to multiple smaller ones lies in being able to do the rough equivalent of for-looping down at a lower language level, which is therefore more efficient than having to use a matlab compiler-generated for-loop structure. I can see no other reason the block diagonal method could win out over simply executing the blocks singly.
> >
> > Roger Stafford (The dreamer)
>
> Yes, you can think of the block diagonal method
> as just avoiding explicit loops, trading them for
> implicit loops in blas calls.
>
> One could also probably have used cellfun or arrayfun
> to solve this problem, but the implicit loops used
> there are not at quite as low a level as the block
> diagonal solution.
>
> John

Subject: overdetermenied system

From: Tim Davis

Date: 15 Jan, 2009 01:10:19

Message: 12 of 12

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gkiegi$alv$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gkiaf5$300$1@fred.mathworks.com>...
> > If your matrix is sparse, then \ really is pretty smart.
> >
> > blocks = rand(4,2,1000);
> > A = blktridiag(blocks,zeros(4,2,999),zeros(4,2,999));
> > rhs = rand(4000,1);
> >
> > tic
> > coef = zeros(2,1000);
> > for i = 1:1000
> > coef(:,i) = blocks(:,:,i)\rhs((1:4) + 4*(i-1));
> > end
> > toc
> > Elapsed time is 0.105407 seconds.
> >
> > tic,coef = A\rhs;toc
> > Elapsed time is 0.044134 seconds.
> >
> > I will admit that the simple loop is pretty fast though,
> > especially if I add in the time to create the sparse
> > block diagonal matrix. blktridiag is not as efficient
> > as it should be here, since it fills in explicit zeros
> > into the upper and lower blocks.
> >
> > John
>
> OK, John and Bruno, I surrender. Your sparse block diagonal method seems to win the day. I was unaware of the 9/14/08 thread mentioned by Bruno.
>
> Nevertheless, it somehow offends my sense of how such an algorithm ought to work. Instead of a procedure that has to be smart enough to ignore off-diagonal zeros in a sparse matrix, I dream of one that never has to deal with them in the first place. Ideally there ought to be a method of carrying out the operations involved in each independent block in a repetitive manner that avoids the inefficiencies of for-loop indexing and its overhead. I suspect the apparent efficiency here of the single backslash operation as opposed to multiple smaller ones lies in being able to do the rough equivalent of for-looping down at a lower language level, which is therefore more efficient than having to use a matlab compiler-generated for-loop structure. I can see no other reason the block diagonal method could win out over simply executing the blocks singly.
>
> Roger Stafford (The dreamer)

There are too many details inside to explain completely ... and I didn't write the sparse QR part of backslash (which is what you're using here).

It should be faster to do a single x=A\b, but backslash does things like "allocate and set to zero an array of size n". That doesn't happen when you solve lots of little blocks for a block diagonal matrix. So I can see cases where backslash does the "wrong thing" because the way the code is written it's assuming it's solving a single system, not lots of teeny ones.

To solve lots of teeny matrices, you really need to use dmperm first (inside backslash). My sparse LU code "KLU" (not UMFPACK) does this.

If you'd like to try a better sparse QR, see SuiteSparseQR in the SuiteSparse package. It might behave differently for this matrix (it's a collection of teeny overdetermined systems, right?).

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