commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Al Chou <>
Subject Re: [math] Does LU decomposition assume matrix is square?
Date Sat, 21 Jun 2003 20:23:34 GMT
--- Phil Steitz <> 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

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

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


Albert Davidson Chou

    Get answers to Mac questions at .

Do you Yahoo!?
SBC Yahoo! DSL - Now only $29.95 per month!

To unsubscribe, e-mail:
For additional commands, e-mail:

View raw message