[math] RealMatrix implementation strategy

[Phil Steitz <ph...@steitz.com>]

I just submitted a patch filling in all but the getRank() methods in 
RealMatrixImpl. I thought a lot about adding a linear solution 
framework, but I decided to implement LU decomposition directly in 
RealMatrixImpl for the following reasons:

1. RealMatrixImpl is an implementation strategy itself.

2. The only thing common among the various strategies for solving linear 
systems is solve(). LU decomposition (what I implemented) supports 
multiple operations (isSingular,getDeterminant,inverse), with the 
decomposition reused once it has been computed.  To implement these 
efficiently in RealMatrixImpl requires a commitment to a solution 
strategy.  I chose LU decompostion using Crout's method with partial 
pivoting (but no scaling).

3. LU Decomp seems to be simply the best general solution.

Other RealMatrix implementations could choose other solution strategies. 
I suppose that we could create an abstract base class including 
implementations of the trivial stuff and then have strategies that 
wanted to use different machinery for solve(), etc., subclass and 
implement.  This may make sense down the road if and when we get other 
RealMatrix implementations.  For now, I would prefer to keep it simple.

There are a couple of numerical issues that I could use some 
confirmation/advice on:

1. How essential is it to add scaling?  I notice that NR does this; but 
other implementations do not.

2. The method that I used to detect singularity was to set a threshold 
for the smallest maximum pivot element during the LU decomposition. 
(this is the element that will end up on the diagonal). I set this, 
somewhat arbitrarily, at 10E-12.  I had it set much lower initially, but 
my singularity tests were failing (singular matrices were coming back as 
nonsingular).  I suppose this should probably be exposed so that the 
user can set it, but I would like to understand first what a good 
default value should be.  Setting it to a high value might avoid 
stability issues with solve() and inverse() for near-singular matrices, 
but it will restrict the domain of application.

Thanks in advance.  I am going out of town for the next week and I may 
have limited or no access to email, so please bear with me if it takes 
me a while to see/respond to feedback.

Phil


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