Thread Subject: Is a QR Decomposition Better than B \ A?

If we have a system of equations Ax=b, then we can solve it using

  x = A \ b

or

  x = pinv(A'*A)*A'*b

or

  [Q,R] = qr(A);
  x = inv(R)*Q'*b;

Which is best?

Also, I thought a QR decomposition always produced a nonsingular R
matrix [meyer, p.313], but Matlab is returning a non-square
matrix. Why?

--Randy

@BOOK{meyer,
  title = "{Matrix Analysis and Applied Linear Algebra}",
  author = "{Carl~D.~Meyer}",
  publisher = "Society for Industrial and Applied Mathematics",
  year = "2000"}

--
% Randy Yates % "The dreamer, the unwoken fool -
%% Fuquay-Varina, NC % in dreams, no pain will kiss the brow..."
%%% 919-577-9882 %
%%%% <yates@ieee.org> % 'Eldorado Overture', *Eldorado*, ELO
http://home.earthlink.net/~yatescr

Randy Yates wrote:
>
>
> If we have a system of equations Ax=b, then we can solve it using
>
> x = A \ b
>
> or
>
> x = pinv(A'*A)*A'*b

NO. NO. NO!!!!

Why do people persist in using this form?

 set('RantMode','off')

Whew. That feels much better now. Forming
the product A'*A is a bad thing to do.
Using pinv does not make it magically
better. Especially when pinv is best
used directly, as in:

 x = pinv(A)*b;

 
> or
>
> [Q,R] = qr(A);
> x = inv(R)*Q'*b;

No, No, No.

 set('RantMode','off')

(Damn that RantMode. It keeps resetting
itself.)

No need to use inv here. It is inefficient
and totally unecessary.

 x = R\*(Q'*b);

Since R is upper triangular, \ is FAR more
efficient here.

 
> Which is best?

Each is different, with different properties.
I talk about some of those differences
in my tips and tricks doc.

 <http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8553&objectType=FILE>

 
> Also, I thought a QR decomposition always produced a nonsingular R
> matrix [meyer, p.313], but Matlab is returning a non-square
> matrix. Why?

No. QR will produce a non-singular
matrix to the extent that your original
matrix is non-singular. That does not
mean that R will be square. It WILL
be upper triangular. Note the difference
between qr(A) and qr(A,0), the economy
version of qr.

A = rand(5,2);
[Q,R] = qr(A)
Q =
     -0.59222 -0.11136 -0.48562 -0.27491 -0.57051
    -0.013877 -0.92778 -0.087834 -0.13611 0.33585
     -0.48139 -0.11727 0.84955 -0.088917 -0.1577
     -0.26814 -0.15247 -0.087584 0.94442 -0.072426
     -0.58775 0.29972 -0.16447 -0.077827 0.72912
R =
      -1.4152 -0.70723
            0 -0.75412
            0 0
            0 0
            0 0

[Q,R] = qr(A,0)
Q =
     -0.59222 -0.11136
    -0.013877 -0.92778
     -0.48139 -0.11727
     -0.26814 -0.15247
     -0.58775 0.29972
R =
      -1.4152 -0.70723
            0 -0.75412

You probably are thinking about the
economy sized version.

HTH,
John

John D'Errico <woodchips@rochester.rr.com> wrote:
> Randy Yates wrote:
>>
>>
>> If we have a system of equations Ax=b, then we can solve it using
>>
[snip]

>
> x = R\*(Q'*b);

A Rantmode induced typo, I think.

x = R \ (Q'*b);

>
> Since R is upper triangular, \ is FAR more
> efficient here.
>
>

Excellent rant!

--
Dr Tristram J. Scott
Energy Consultant

"John D'Errico" <woodchips@rochester.rr.com> writes:

> Randy Yates wrote:
>>
>>
>> If we have a system of equations Ax=b, then we can solve it using
>>
>> x = A \ b
>>
>> or
>>
>> x = pinv(A'*A)*A'*b
>
> NO. NO. NO!!!!
>
> Why do people persist in using this form?

Because it's what you see in the linear algebra books and it's so cool
in theory?

To Meyer's credit, he mentions that it's a bad thing to do numerically
in his text and had the QR method as a suggested alternative. Indeed
that's the reason I'm asking here now.

I think his book is super excellent, by the way.

> set('RantMode','off')
>
> Whew. That feels much better now. Forming
> the product A'*A is a bad thing to do.
> Using pinv does not make it magically
> better. Especially when pinv is best
> used directly, as in:
>
> x = pinv(A)*b;

Hmmm. Okie-dokie. Sorry,

>> or
>>
>> [Q,R] = qr(A);
>> x = inv(R)*Q'*b;
>
> No, No, No.
>
> set('RantMode','off')
>
> (Damn that RantMode. It keeps resetting
> itself.)

Yeah, damn-it.

> No need to use inv here. It is inefficient
> and totally unecessary.
>
> x = R\*(Q'*b);
>
> Since R is upper triangular, \ is FAR more
> efficient here.

Not only that, it runs a lot faster... :)

>> Which is best?
>
> Each is different, with different properties.
> I talk about some of those differences
> in my tips and tricks doc.
>
> <http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8553&objectType=FILE>

I appreciate the resource, but I must admit the format is so
unconventional that I had a hard time viewing it even using
Matlab. Why not prepare a simple pdf instead of embedding
documentation (including specific font information) in the m
files?

>> Also, I thought a QR decomposition always produced a nonsingular R
>> matrix [meyer, p.313], but Matlab is returning a non-square
>> matrix. Why?
>
> No. QR will produce a non-singular
> matrix to the extent that your original
> matrix is non-singular. That does not
> mean that R will be square. It WILL
> be upper triangular. Note the difference
> between qr(A) and qr(A,0), the economy
> version of qr.

OK, got it - thanks.
--
% Randy Yates % "Though you ride on the wheels of tomorrow,
%% Fuquay-Varina, NC % you still wander the fields of your
%%% 919-577-9882 % sorrow."
%%%% <yates@ieee.org> % '21st Century Man', *Time*, ELO
http://home.earthlink.net/~yatescr

If we have the equation Ax=b.
A is sparse.
The number of x is smaller than the number of equation.
Is x=A\b the best and quickest way to solve x ?

Randy Yates wrote:
>
>
> "John D'Errico" <woodchips@rochester.rr.com> writes:
>
>> Randy Yates wrote:
>>>
>>>
>>> If we have a system of equations Ax=b, then we can solve it
> using
>>>
>>> x = A \ b
>>>
>>> or
>>>
>>> x = pinv(A'*A)*A'*b
>>
>> NO. NO. NO!!!!
>>
>> Why do people persist in using this form?
>
> Because it's what you see in the linear algebra books and it's so
> cool
> in theory?

No, actually its pretty poor in theory,
once you get into the theory.

And, I seriously hope that you will
not find an expression like

  pinv(A'*A)*A'*b

in any recent linear algebra text
worth reading. This shows a lack of
understanding of what pinv does. Why
is it poor? Pinv uses a SVD. But if
you were going to go to the effort
of computing a SVD, applying it to
A'*A squares the singular values.
This loses lots of precision, half
your digits go away into the bit
bucket. At this point there seems
little purpose in throwing a
pseudo-inverse at it. If you were
worried about singularity before,
then why square the singular values?

Yes, old texts do tend to have the
normal equations form:

  inv(A'*A)*A'*b

As I said, its a poor thing to do
with inv or pinv.

>
> To Meyer's credit, he mentions that it's a bad thing to do
> numerically
> in his text and had the QR method as a suggested alternative.

Good.

(snip)

>>> Which is best?
>>
>> Each is different, with different properties.
>> I talk about some of those differences
>> in my tips and tricks doc.
>>
>> <http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8553&objectType=FILE>
>
> I appreciate the resource, but I must admit the format is so
> unconventional that I had a hard time viewing it even using
> Matlab. Why not prepare a simple pdf instead of embedding
> documentation (including specific font information) in the m
> files?

The idea is as an m-file you can change
things and see what happens. Its fully
executable in Matlab. I've considered
moving it to LaTeX/pdf, but then its
just another book to read.

John

"John D'Errico" <woodchips@rochester.rr.com> writes:
> [...]

Hey, dude: you're an arrogant asshole. Your arrogance offsets any
desire I might have had to learn from you.

Try coming down a few notches if you really want to help folks.
--
% Randy Yates % "...the answer lies within your soul
%% Fuquay-Varina, NC % 'cause no one knows which side
%%% 919-577-9882 % the coin will fall."
%%%% <yates@ieee.org> % 'Big Wheels', *Out of the Blue*, ELO
http://home.earthlink.net/~yatescr

Randy Yates wrote:
> If we have a system of equations Ax=b, then we can solve it using
>
> x = A \ b
>
> or
>
> x = pinv(A'*A)*A'*b
>
> or
>
> [Q,R] = qr(A);
> x = inv(R)*Q'*b;
>
> Which is best?
>
> Also, I thought a QR decomposition always produced a nonsingular R
> matrix [meyer, p.313], but Matlab is returning a non-square
> matrix. Why?
>
> --Randy
>
> @BOOK{meyer,
> title = "{Matrix Analysis and Applied Linear Algebra}",
> author = "{Carl~D.~Meyer}",
> publisher = "Society for Industrial and Applied Mathematics",
> year = "2000"}
>

Which alternative is "best" will likely depend on your application. The
matlab function reference pages for mldivide, pinv, and qr will give you
excruciating detail about the differences (and similarities) between the
three. In a nutshell, backslash-division will try to intelligently
select a solution algorithm based on e.g. whether A is square or not,
and sparse or not (and for non-square matrices it uses the QR
decomposition) and will produce a solution with the minimum number of
non-zero entries. Using pinv (in the form x = pinv(A) * b; the
equation you give works too but why bother with the extra operations?)
will give a solution with minimum norm. Unless you know you want this
I'd suggest using \ and letting Matlab do the thinking about which
algorithm to use.

As to your second question, it is the Q matrix which is guaranteed to be
non-singular (even better, orthogonal). R must always have the same
size as A, and be singular if and only if A is.

Robert E. Beaudoin

"Robert E. Beaudoin" <rbeaudoin@comcast.net> writes:

> Randy Yates wrote:
>> If we have a system of equations Ax=b, then we can solve it using
>> x = A \ b or x = pinv(A'*A)*A'*b
>> or
>> [Q,R] = qr(A);
>> x = inv(R)*Q'*b;
>> Which is best? Also, I thought a QR decomposition always produced a
>> nonsingular R
>> matrix [meyer, p.313], but Matlab is returning a non-square
>> matrix. Why?
>> --Randy
>> @BOOK{meyer,
>> title = "{Matrix Analysis and Applied Linear Algebra}",
>> author = "{Carl~D.~Meyer}",
>> publisher = "Society for Industrial and Applied Mathematics",
>> year = "2000"}
>>
>
> Which alternative is "best" will likely depend on your application.
> The matlab function reference pages for mldivide, pinv, and qr will
> give you excruciating detail about the differences (and similarities)
> between the three. In a nutshell, backslash-division will try to
> intelligently select a solution algorithm based on e.g. whether A is
> square or not, and sparse or not (and for non-square matrices it uses
> the QR decomposition) and will produce a solution with the minimum
> number of non-zero entries. Using pinv (in the form x = pinv(A) * b;
> the equation you give works too but why bother with the extra
> operations?) will give a solution with minimum norm. Unless you know
> you want this I'd suggest using \ and letting Matlab do the thinking
> about which algorithm to use.
>
> As to your second question, it is the Q matrix which is guaranteed to
> be non-singular (even better, orthogonal). R must always have the
> same size as A, and be singular if and only if A is.

Hello Robert,

Thank you for your help. I took another look at the matlab help page
for qr and I don't see any discussion of the differences or similarities.
the mldivide page does have some information. Perhaps I'm using an older
version than you (6.1.0.450 (R12.1))?

What I'm trying to ask is, which method gives me the minimum
squared error (minimum norm of residual)?
--
% Randy Yates % "Remember the good old 1980's, when
%% Fuquay-Varina, NC % things were so uncomplicated?"
%%% 919-577-9882 % 'Ticket To The Moon'
%%%% <yates@ieee.org> % *Time*, Electric Light Orchestra
http://home.earthlink.net/~yatescr

Wow. Randy.

Reread what John typed - I think you're reading in something that he didn't put in.

Scott

On Jul 16, 8:38 am, Randy Yates <ya...@ieee.org> wrote:
> "Robert E. Beaudoin" <rbeaud...@comcast.net> writes:
>
>
>
>
>
> > Randy Yates wrote:
> >> If we have a system of equations Ax=b, then we can solve it using
> >> x = A \ b or x = pinv(A'*A)*A'*b
> >> or
> >> [Q,R] = qr(A);
> >> x = inv(R)*Q'*b;
> >> Which is best? Also, I thought a QR decomposition always produced a
> >> nonsingular R
> >> matrix [meyer, p.313], but Matlab is returning a non-square
> >> matrix. Why?
> >> --Randy
> >> @BOOK{meyer,
> >> title = "{Matrix Analysis and Applied Linear Algebra}",
> >> author = "{Carl~D.~Meyer}",
> >> publisher = "Society for Industrial and Applied Mathematics",
> >> year = "2000"}
>
> > Which alternative is "best" will likely depend on your application.
> > The matlab function reference pages for mldivide, pinv, and qr will
> > give you excruciating detail about the differences (and similarities)
> > between the three. In a nutshell, backslash-division will try to
> > intelligently select a solution algorithm based on e.g. whether A is
> > square or not, and sparse or not (and for non-square matrices it uses
> > the QR decomposition) and will produce a solution with the minimum
> > number of non-zero entries. Using pinv (in the form x = pinv(A) * b;
> > the equation you give works too but why bother with the extra
> > operations?) will give a solution with minimum norm. Unless you know
> > you want this I'd suggest using \ and letting Matlab do the thinking
> > about which algorithm to use.
>
> > As to your second question, it is the Q matrix which is guaranteed to
> > be non-singular (even better, orthogonal). R must always have the
> > same size as A, and be singular if and only if A is.
>
> Hello Robert,
>
> Thank you for your help. I took another look at the matlab help page
> for qr and I don't see any discussion of the differences or similarities.
> the mldivide page does have some information. Perhaps I'm using an older
> version than you (6.1.0.450 (R12.1))?
>
> What I'm trying to ask is, which method gives me the minimum
> squared error (minimum norm of residual)?
> --
> % Randy Yates % "Remember the good old 1980's, when
> %% Fuquay-Varina, NC % things were so uncomplicated?"
> %%% 919-577-9882 % 'Ticket To The Moon'
> %%%% <ya...@ieee.org> % *Time*, Electric Light Orchestrahttp://home.earthlink.net/~yatescr- Hide quoted text -
>
> - Show quoted text -

Hello Randy,

Are you looking for the least squared error solution to Ax=b, then
solve

x = (((A^t)A)^-1) *b

where A^t is A transpose and ()^-1 is inverse

if you let A = QR then you find


x = (R^-1)(Q^t)b

I hope this helps

Clay

On Jul 16, 3:54 pm, Clay <phys...@bellsouth.net> wrote:
> On Jul 16, 8:38 am, Randy Yates <ya...@ieee.org> wrote:
>
>
>
>
>
> > "Robert E. Beaudoin" <rbeaud...@comcast.net> writes:
>
> > > Randy Yates wrote:
> > >> If we have a system of equations Ax=b, then we can solve it using
> > >> x = A \ b or x = pinv(A'*A)*A'*b
> > >> or
> > >> [Q,R] = qr(A);
> > >> x = inv(R)*Q'*b;
> > >> Which is best? Also, I thought a QR decomposition always produced a
> > >> nonsingular R
> > >> matrix [meyer, p.313], but Matlab is returning a non-square
> > >> matrix. Why?
> > >> --Randy
> > >> @BOOK{meyer,
> > >> title = "{Matrix Analysis and Applied Linear Algebra}",
> > >> author = "{Carl~D.~Meyer}",
> > >> publisher = "Society for Industrial and Applied Mathematics",
> > >> year = "2000"}
>
> > > Which alternative is "best" will likely depend on your application.
> > > The matlab function reference pages for mldivide, pinv, and qr will
> > > give you excruciating detail about the differences (and similarities)
> > > between the three. In a nutshell, backslash-division will try to
> > > intelligently select a solution algorithm based on e.g. whether A is
> > > square or not, and sparse or not (and for non-square matrices it uses
> > > the QR decomposition) and will produce a solution with the minimum
> > > number of non-zero entries. Using pinv (in the form x = pinv(A) * b;
> > > the equation you give works too but why bother with the extra
> > > operations?) will give a solution with minimum norm. Unless you know
> > > you want this I'd suggest using \ and letting Matlab do the thinking
> > > about which algorithm to use.
>
> > > As to your second question, it is the Q matrix which is guaranteed to
> > > be non-singular (even better, orthogonal). R must always have the
> > > same size as A, and be singular if and only if A is.
>
> > Hello Robert,
>
> > Thank you for your help. I took another look at the matlab help page
> > for qr and I don't see any discussion of the differences or similarities.
> > the mldivide page does have some information. Perhaps I'm using an older
> > version than you (6.1.0.450 (R12.1))?
>
> > What I'm trying to ask is, which method gives me the minimum
> > squared error (minimum norm of residual)?
> > --
> > % Randy Yates % "Remember the good old 1980's, when
> > %% Fuquay-Varina, NC % things were so uncomplicated?"
> > %%% 919-577-9882 % 'Ticket To The Moon'
> > %%%% <ya...@ieee.org> % *Time*, Electric Light Orchestrahttp://home.earthlink.net/~yatescr-Hide quoted text -
>
> > - Show quoted text -
>
> Hello Randy,
>
> Are you looking for the least squared error solution to Ax=b, then
> solve
>
> x = (((A^t)A)^-1) *b
>
> where A^t is A transpose and ()^-1 is inverse
>
> if you let A = QR then you find
>
> x = (R^-1)(Q^t)b
>
> I hope this helps
>
> Clay- Hide quoted text -
>
> - Show quoted text -


Sorry I didn't check it fully before sending

x = (((A^t)A)^-1)(A^t)b

or

x = (R^-1)(Q^t)b

will give LSQR solutions to Ax=b

Clay

Clay <physics@bellsouth.net> writes:
> [...]
> Sorry I didn't check it fully before sending
>
> x = (((A^t)A)^-1)(A^t)b
>
> or
>
> x = (R^-1)(Q^t)b
>
> will give LSQR solutions to Ax=b
>
> Clay

As Robert said, that's only good if A is non-singular, and it ain't
necessarily in my application.

I know these are all least-squares solutions. I'm asking which ones are
"better." That is, which one provides a smaller squared-error?
--
% Randy Yates % "She's sweet on Wagner-I think she'd die for Beethoven.
%% Fuquay-Varina, NC % She love the way Puccini lays down a tune, and
%%% 919-577-9882 % Verdi's always creepin' from her room."
%%%% <yates@ieee.org> % "Rockaria", *A New World Record*, ELO
http://home.earthlink.net/~yatescr

Randy Yates wrote:
> Clay <physics@bellsouth.net> writes:
>> [...]
>> Sorry I didn't check it fully before sending
>>
>> x = (((A^t)A)^-1)(A^t)b
>>
>> or
>>
>> x = (R^-1)(Q^t)b
>>
>> will give LSQR solutions to Ax=b
>>
>> Clay
>
> As Robert said, that's only good if A is non-singular, and it ain't
> necessarily in my application.
>
> I know these are all least-squares solutions. I'm asking which ones are
> "better." That is, which one provides a smaller squared-error?

According to the Matlab help pages (I'm looking at them via the web as I
  don't have a Matlab installation at home, so I think they apply to the
latest Matlab release) either mldivide or pinv will return a
least-square solution. That means they both return a solution
minimizing the (squared) error. So when the least-square solution is
unique they will return the same answer (within round-off error). But
if A is n-by-m with n > m > rank(A) there can be more than one solution
minimizing the squared error; in that case A\b returns a least-square
solution with no more nonzero entries than the rank of A, whereas
pinv(A) * b is a least-square solution which itself has minimum norm
(among all least-square solutions). The Mathworks docs give the example

A = [ 1 2 3
       4 5 6
       7 8 9
      10 11 12]

b = [2 5 8 11]'.

The minimum squared error is zero, an obvious least-square (i.e. exact)
solution is [0 1 0]', but A\b yields the least-square solution
[1/2 0 1/2]' whereas pinv(A) * b yields the least-square solution
[1/3 1/3 1/3]'.

Normally you wouldn't really care which of these you get; as a practical
matter A\b is likely to be faster than pinv(A) * b but could give you
annoying warning messages if A is rank-deficient (as in this example).
Depending on your application this rank-deficient case may never arise.

Hope that helps. I'll be traveling for the rest of this week, so I
won't be able to answer any further questions for a while.

Robert E. Beaudoin

"Robert E. Beaudoin" <rbeaudoin@comcast.net> writes:
> [...]
> Hope that helps. I'll be traveling for the rest of this week, so I
> won't be able to answer any further questions for a while.

Yes, that helps immensely. Thank you, Robert.
--
% Randy Yates % "...the answer lies within your soul
%% Fuquay-Varina, NC % 'cause no one knows which side
%%% 919-577-9882 % the coin will fall."
%%%% <yates@ieee.org> % 'Big Wheels', *Out of the Blue*, ELO
http://home.earthlink.net/~yatescr

Dear All,

FYI,

The following is a copy from the linear algebra course of Prof.
Gilbert Strang at MIT on backslash command to solve Ax=b. (http://www-
math.mit.edu/18085/). The pdf file name is "backslash.pdf"

Hope this helps.

Sangdon Lee, Ph.D.
GM Tech. Center

%%%%% start %%%%%%%%%%

MATLAB's backslash command to solve Ax = b

 x=A\b for dense A performs these steps (stopping when successful):

1. If A is upper or lower triangular, solve by back/forward
substitution
2. If A is permutation of triangular matrix, solve by permuted back
substitution
(useful for [L,U]=lu(A) since L is permuted)
3. If A is symmetric/hermitian
- Check if all diagonal elements are positive
- Try Cholesky, if successful solve by back substitutions
4. If A is Hessenberg (upper triangular plus one subdiagonal), reduce
to upper tri-
angular then solve by back substitution
5. If A is square, factorize PA = LU and solve by back substitutions
6. If A is not square, run Householder QR, solve least squares
problem.

Mathworks documentation:
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/mldivide.html#1002049

%%%% end %%%%%%

Randy Yates wrote:
>
>
> "John D'Errico" <woodchips@rochester.rr.com> writes:
>> [...]
>
> Hey, dude: you're an arrogant asshole. Your arrogance offsets any
> desire I might have had to learn from you.
>
> Try coming down a few notches if you really want to help folks.
> --
> % Randy Yates % "...the answer lies within your
> soul
> %% Fuquay-Varina, NC % 'cause no one knows which
> side
> %%% 919-577-9882 % the coin will
> fall."
> %%%% <yates@ieee.org> % 'Big Wheels', *Out of the
Blue*,
> ELO
> <http://home.earthlink.net/~yatescr>
>

I'd like to thank John for all his contributions to the forum and for
tips and tricks etc. I can tell from the download record, other
apprecriate his work as well. I, for one, learn a lot reading his
posts.

"Randy Yates" <yates@ieee.org> wrote in message
news:m3abtz9895.fsf@ieee.org...
> "John D'Errico" <woodchips@rochester.rr.com> writes:
>> [...]
>
> Hey, dude: you're an arrogant asshole. Your arrogance offsets any
> desire I might have had to learn from you.

Great... then please go away!

By the way, I saw nothing arrogant but much that I read as knowledgable.

> Try coming down a few notches if you really want to help folks.
> --
> % Randy Yates % "...the answer lies within your soul
> %% Fuquay-Varina, NC % 'cause no one knows which side
> %%% 919-577-9882 % the coin will fall."
> %%%% <yates@ieee.org> % 'Big Wheels', *Out of the Blue*, ELO
> http://home.earthlink.net/~yatescr