You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by Al Chou <ho...@yahoo.com> on 2003/06/21 19:55:06 UTC

[math] Does LU decomposition assume matrix is square?

I finally decided that cubic spline would be my first attempt at implementing
interpolation, partly because of the difficulty of finding an alternative
reference to NR for the Stoer & Bulirsch rational function method, partly
because I started to have doubts about the desirability of interpolation using
high-order polynomials/rationals, and partly because I think being forced to
implement tridiagonal linear systems solution is probably good for the library.
 I could go the NR route and embed the tridiagonal solver into the cubic spline
routine, but I think it's worthwhile to provide the tridiagonal solver
separately for others to use, given the frequency with which tridiagonal
systems appear (in physics anyway; am I dreaming too much to think that someone
might use this library to solve second order differential equations by finite
differences?).

I happened to notice as I started to plan my work that RealMatrixImpl.solve,
which uses LU decomposition, doesn't explicitly check whether the matrix is
square before proceeding with the decomposition.  I think (but am not sure)
that LU decomposition assumes the matrix is square.  Currently it already
checks whether the vector of right-hand-sides of the equations is equal to the
number of rows in the matrix, which I think implicitly assumes that the matrix
is square (otherwise it would probably be more correct to check against the
number of columns).


Al

=====
Albert Davidson Chou

    Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
SBC Yahoo! DSL - Now only $29.95 per month!
http://sbc.yahoo.com

---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Re: [math] Does LU decomposition assume matrix is square?

Posted by Al Chou <ho...@yahoo.com>.
--- Phil Steitz <ph...@steitz.com> wrote:
> Al Chou wrote:
> > I finally decided that cubic spline would be my first attempt at
> implementing
> > interpolation, partly because of the difficulty of finding an alternative
> > reference to NR for the Stoer & Bulirsch rational function method, partly
> > because I started to have doubts about the desirability of interpolation
> using
> > high-order polynomials/rationals
> 
> +1  This is also easier to document and understand for users.

Yes, in reacquainting myself with the algorithms I rediscovered one of the
subtleties of using polynomial or rational interpolation, which is that you
probably should pick a small subset of your tabulated points near where you
want to interpolate and use only those in the interpolation, thus keeping the
degree of the interpolating function kind of low and hopefully limiting the
arbitrary wiggliness of the interpolating function.  Cubic spline allows you to
do what seems like the easiest thing, which is provide a whole table of data
points and ask for interpolation anywhere within the tabulated domain.  Of
course, you can do a similar thing with polynomial or rational interpolation,
but either the library or the user would have to provide the wrapper that asks
for how many points to use in the neighborhood of the interpolation point, or
alternatively the degree of the interpolating function to use, both of which
require more decisions on the user's part in the process of using
interpolation.


> , and partly because I think being forced to
> > implement tridiagonal linear systems solution is probably good for the
> library.
> > I could go the NR route and embed the tridiagonal solver into the cubic
> spline
> > routine, but I think it's worthwhile to provide the tridiagonal solver
> > separately for others to use, given the frequency with which tridiagonal
> > systems appear (in physics anyway; am I dreaming too much to think that
> someone
> > might use this library to solve second order differential equations by
> finite
> > differences?).
> 
> I would either embed the tridiagonal solution or just use
> RealMatrixImpl.solve(). Personally, I would probably take the second 
> approach, for initial release, since the implementation exists.  If 
> users complain about the speed (which will only happen if they are 
> fitting splines using large numbers of points), we can modify the 
> implementation to use the faster approach. I would not see a tridiagonal 
> solver as something that we should add to intial release of commons-math.

Good suggestion.  I'll start from there, given that I discovered more design
decisions to be made as I sketched converting RealMatrixImpl to a
RealTridiagonalMatrix....


> > I happened to notice as I started to plan my work that
> RealMatrixImpl.solve,
> > which uses LU decomposition, doesn't explicitly check whether the matrix is
> > square before proceeding with the decomposition.  I think (but am not sure)
> > that LU decomposition assumes the matrix is square.
> 
> No, but it does require that the number of rows >= the number of 
> columns.  I will add this check.  Good catch.
> 
> >  Currently it already
> > checks whether the vector of right-hand-sides of the equations is equal to
> the
> > number of rows in the matrix, which I think implicitly assumes that the
> matrix
> > is square (otherwise it would probably be more correct to check against the
> > number of columns).
> 
> The check against the constant matrix column dimension is to make sure 
> that the linear system(s) represented by the matrix equation are 
> well-defined.  You are correct that the coefficient matrix also needs to 
> be square. I will add this to solve().  Once again, good catch!

Cool, glad I haven't lost (all) my marbles. <g>



Al

=====
Albert Davidson Chou

    Get answers to Mac questions at http://www.Mac-Mgrs.org/ .

__________________________________
Do you Yahoo!?
SBC Yahoo! DSL - Now only $29.95 per month!
http://sbc.yahoo.com

---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Re: [math] Does LU decomposition assume matrix is square?

Posted by Phil Steitz <ph...@steitz.com>.
Al Chou wrote:
> I finally decided that cubic spline would be my first attempt at implementing
> interpolation, partly because of the difficulty of finding an alternative
> reference to NR for the Stoer & Bulirsch rational function method, partly
> because I started to have doubts about the desirability of interpolation using
> high-order polynomials/rationals

+1  This is also easier to document and understand for users.

, and partly because I think being forced to
> implement tridiagonal linear systems solution is probably good for the library.
> I could go the NR route and embed the tridiagonal solver into the cubic spline
> routine, but I think it's worthwhile to provide the tridiagonal solver
> separately for others to use, given the frequency with which tridiagonal
> systems appear (in physics anyway; am I dreaming too much to think that someone
> might use this library to solve second order differential equations by finite
> differences?).

I would either embed the tridiagonal solution or just use
RealMatrixImpl.solve(). Personally, I would probably take the second 
approach, for initial release, since the implementation exists.  If 
users complain about the speed (which will only happen if they are 
fitting splines using large numbers of points), we can modify the 
implementation to use the faster approach. I would not see a tridiagonal 
solver as something that we should add to intial release of commons-math.

> 
> I happened to notice as I started to plan my work that RealMatrixImpl.solve,
> which uses LU decomposition, doesn't explicitly check whether the matrix is
> square before proceeding with the decomposition.  I think (but am not sure)
> that LU decomposition assumes the matrix is square.

No, but it does require that the number of rows >= the number of 
columns.  I will add this check.  Good catch.

>  Currently it already
> checks whether the vector of right-hand-sides of the equations is equal to the
> number of rows in the matrix, which I think implicitly assumes that the matrix
> is square (otherwise it would probably be more correct to check against the
> number of columns).

The check against the constant matrix column dimension is to make sure 
that the linear system(s) represented by the matrix equation are 
well-defined.  You are correct that the coefficient matrix also needs to 
be square. I will add this to solve().  Once again, good catch!

> 
> 
> Al
> 
> =====
> Albert Davidson Chou
> 
>     Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
> 
> __________________________________
> Do you Yahoo!?
> SBC Yahoo! DSL - Now only $29.95 per month!
> http://sbc.yahoo.com
> 
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
> For additional commands, e-mail: commons-dev-help@jakarta.apache.org
> 




---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org