I have an initial implementation of LSMR that I will be checking into the
Mahout project shortly. If you want to have such an algorithm in commons
math, you would be entirely free to grab that and adapt it for possible
inclusion in commons math (or your own use).
2011/4/12 Sébastien Brisard <sebastien.brisard@m4x.org>
> Dear All,
> first of all, I should say that I'm new to the CommonsMath community, and
> I
> therefore apologize for any mistake/inconsistency.
> I'm currently developing a numerical method for the simulation of highly
> heterogeneous materials. This requires the inversion of a (big) linear
> system,
> whose matrix A is *not* known in closed form. The only thing I know about A
> is
> how to compute the product A.x for any vector x. True enough, this would
> allow
> me to compute the coefficients of the matrix, but
> 1. I do not need them,
> 2. It would waste memory space.
>
> Many iterative solvers actually do not require the coefficients of matrix
> A.
> For example, I've implemented CG and SYMMLQ, only assuming that A is
> defined
> by a method operate(double[] arg, double[] img).
>
> If you feel that this code might find its place in CommonsMath, I'm very
> willing to contribute it. However, here are a few questions
> 1. The matrix as I defined it does not fit into the
> org.apache.commons.math.linear.RealMatrix interface (since I do not want to
> define add, getRow, etc...). For the time being, I defined a so called
> LinearMap interface, which inherits from AnyMatrix, but only defines one
> additional method
> void operate(double[] arg, double[] img)
>
> 2. BTW, the vector in which the image A.x of the argument x must be stored
> is
> passed as an argument, in order to limit memory allocation. I always do
> that,
> but noticed that it is not commonpractice in CommonsMath. Should I give
> up
> this habit? Is it badpractice? Is the garbage collector so efficient
> nowadays
> that such precaution has become pointless?
>
> 2bis. If you think I should substitute
> double[] operate(double[] arg)
> to
> void operate(double[] arg, double[] img)
> then RealMatrix might extend LinearMap (or whatever its name).
>
> 3. Regarding preconditionners: the GC or SYMMLQ solvers can take as an
> argument a preconditionner M. Again, we do not need to know its
> coefficients.
> The only thing which is required is to compute M.x or inv(M).y. So I have
> defined another interface, namely InvertibleLinearMap, which inherits from
> LinearMap, with one additional method
> void solve(double[] b, double[] x)
> Any thoughts?
>
> 3. Regarding potential exceptions that may be returned. I'm puzzled by the
> constructors of the exceptions which might be of use for me, namely
> NonSymmetricMatrixException
> NonPositiveDefiniteMatrixException
>
> both exceptions refer to row/column indices, which I do not know! For
> example,
> I would like to return a NonPositiveDefiniteMatrixException if, during the
> iterations, I find that transpose(x).A.x < 0, without referring to a
> specific
> row or column. Besides, such an exception might be thrown either when M or
> A
> is nonsymmetric (both passed to the GC/SYMMLQ solver at the same time),
> and I
> would like the message of the exception to state wether the problem came
> from
> the matrix A itself, or the preconditionner M. How should I do? Should I
> define my own exceptions? This seems to me to be a waste...
>
> I'm looking forward to having the opportunity to sharing my code.
> Meanwhile, I
> do apologize if the previous remarks and questions sound silly.
>
> Best,
> Sebastien
>
> 
> To unsubscribe, email: devunsubscribe@commons.apache.org
> For additional commands, email: devhelp@commons.apache.org
>
>
