Prev: Contour maps
Next: Integration of Simulink Model Compare in SVN
From: Rune Allnor on 4 Oct 2009 09:21
On 4 Okt, 02:45, "leo nidas" <bleonida...(a)yahoo.gr> wrote:
> Hi there,
>
> I want to solve or b the : (X'X)*b=X'*y. I have some issues with the inverse of X'*X.

Find a book on elementary linear algebra, read it,
contemplate it, and then use the ages-old standard way
to handle these things:

b = inv(X'X)X'y

Compute the SVD of X:

X = USV where S = diagonal, U'U = UU' =I, V'V = VV'= I

Insert in the expression to find

b = V'inv(S)U'y.

Since S is diagonal, S^-1 is easily computed.

Now, the matlab function PINV does all this,
so if you want a canned routine you can use it.
However, by writing this out in full, you can get
a useful solution by inspecting the elements on
the diagonal of S:

If |s(n,n)| ~ 0 for some n, the corresponding
term does not contribute in the data model.
However, the corresponding term 1/s(n,n) totally
dominates the pseudo inverse.

So in the pseudo inverse you will want to skip
all terms with singular values smaller than some
tolerance.

PINV does this, with the TOL argument that you
as user have to provide in advance - I suspect this
is an absolute tolerance. I would prefer a relative
tolerance, like

RTOL = eps('double')*max(diag(S));

in which case you will have to roll your own.

Rune
From: Petros Mpogiatzis on 4 Oct 2009 09:45
"Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <haa48p$in7$1(a)fred.mathworks.com>...
> "Tim Davis" <davis(a)cise.ufl.edu> wrote in message <haa339$1tk$1(a)fred.mathworks.com>...
>
> > And you should never (well, almost never) do x=(A'*A)\(A'*b) to solve a least squares problem.
>
> Right, but what you could do is
>
> x = (A'*A + L'*L) \ (A'*b)
>
> where L is appropriate Tikhonov matrix. In the standard form L could be sqrt(lambda)*speye(n).
>
> It is anyway dangerous not using any kind of regularization while solving system involving a ill-conditioned matrix (it is *not* the case of OP's problem, since the problem he had is simply due to a bad scaling).
>
> Bruno

or redefine your initial system to include the L'L constraint and then to solve the redefined system
x=Aext\bext
From: Bruno Luong on 4 Oct 2009 10:05
Rune Allnor <allnor(a)tele.ntnu.no> wrote in message <6ed887e5-8dc1-4d0d-8d96-ef4ac3c7b5c7(a)j19g2000vbp.googlegroups.com>...
> On 4 Okt, 02:45, "leo nidas" <bleonida...(a)yahoo.gr> wrote:

>
> PINV does this, with the TOL argument that you
> as user have to provide in advance - I suspect this
> is an absolute tolerance. I would prefer a relative
> tolerance, like
>

Using PINV has two drawbacks:
1. It call for an expansive SVD
2. It does not work on sparse matrix

That's the reason I wrote a PSEUDOINVERSE on FEX based on QR factorization. Admittedly the factorization will be an expensive strategy for sparse matrix in term of storing space, and it is probably is impossible to apply on large scale problems (such as dimensions > 10e4, and yes there is a TONE of applications that need to solve such large system - forward or backward).

I might complete with an iterative regularized solver soon, even if it's no longer a "pseudo inverse" however both the practical scope of both have large overlapping.

> PINV does this, with the TOL argument that you
> as user have to provide in advance - I suspect this
> is an absolute tolerance.

Correct.

>I would prefer a relative
> tolerance, like
>
> RTOL = eps('double')*max(diag(S));

No big deal, just pass input TOL as the desired relative tolerance multiplied by norm/normest of the input matrix (which is the largest s(i,i)).

Bruno
From: Rune Allnor on 4 Oct 2009 10:24
On 4 Okt, 16:05, "Bruno Luong" <b.lu...(a)fogale.findmycountry> wrote:
> Rune Allnor <all...(a)tele.ntnu.no> wrote in message <6ed887e5-8dc1-4d0d-8d96-ef4ac3c7b...(a)j19g2000vbp.googlegroups.com>...
> > On 4 Okt, 02:45, "leo nidas" <bleonida...(a)yahoo.gr> wrote:
>
> > PINV does this, with the TOL argument that you
> > as user have to provide in advance - I suspect this
> > is an absolute tolerance. I would prefer a relative
> > tolerance, like
>
> Using PINV has two drawbacks:
> 1. It call for an expansive SVD

First come up with something that actually works,
only then criticize speed.

> 2. It does not work on sparse matrix

Were there any sparse matrices in this thread?

> > PINV does this, with the TOL argument that you
> > as user have to provide in advance - I suspect this
> > is an absolute tolerance.
>
> Correct.
>
> >I would prefer a relative
> > tolerance, like
>
> > RTOL = eps('double')*max(diag(S));
>
> No big deal, just pass input TOL as the desired relative tolerance multiplied by norm/normest of the input matrix (which is the largest s(i,i)).

....which requires an equally expensive SVD to compute the
norm before one calls the SVD to compute the PINV...

Get a grip!

Rune

From: Bruno Luong on 4 Oct 2009 10:46
Rune Allnor <allnor(a)tele.ntnu.no> wrote in message <0d69cf8a-2563-4458-b112-401524512d0d(a)p36g2000vbn.googlegroups.com>...

>
> ...which requires an equally expensive SVD to compute the
> norm before one calls the SVD to compute the PINV...
>

WRONG again, the norm estimation use Krylov space technique (Arnoldi iteration). With very few iterations on Rayleigh ratio will converge toward the norm. There is a difference between finding FEW eigenvalue (here once) and all of them.

Your turn to get a grip and open your old algebra book Rune.

Bruno
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4
Prev: Contour maps
Next: Integration of Simulink Model Compare in SVN