Re: [eigen] solve API

[Manoj Rajagopalan <rmanoj@xxxxxxxxx>]

On Monday 28 June 2010 05:42:49 am Gael Guennebaud wrote:

> Seriously, currently our solvers only allows you to compute x = A^1 b
> via some decompositions of A. However, it is sometime needed to also
> solve for x = A^-T b or x = b A^1.

==== Great! This is very useful in building (left-)right-preconditioned sparse solvers using Eigen as
the implementation.

Here are my 2 bits. I am still learning how Eigen3 is implemented so some of my suggestions
might not be implementable because of difficulties with template expressions but please educate me.


[I] How about free functions: (with appropriate template params)

   MatrixBase solve(MatrixBase const& A, MatrixBase const& b); // on-left
   MatrixBase solve(MatrixBase const& b, MatrixBase const& A); // on-right

  Usage:

    x = solve(A,b); // left
    x = solve(b,A); // right
    x = solve(b, A.triangularView<Lower>()); // Gael's example

** Are such "free" functions what you call imperative API in eigen-speak?

**The solve() function would return some matrix-expression (a "proxy" in eigen-speak?)
that would get evaluated and compiled into the appropriate left- or right- solution when
the = expression is evaluated. At that time pointer checks would resolve if y and x are
aliased so that in-place solve could be called instead?

** The issue with this approach is which method to use for solving? LU or other? With A.lu().solve()
    this issue is pre-empted.

** We could rename the free solve() function to lu_solve() to circumvent this issue. We could also provide
    an optional template parameter to specify full or partial pivoting.

** This would still, internally, use the Eigen class API to implement the actual solution (Gael's/Benoit's suggestions)
    but the user would remain insulated from any such changes. The free functions would export a standard,
    readable interface and hide the details from the user.

** For "b" being a complex expression we could have a free solve_in_place()
    // on left
    solve_in_place(A.selfadjointView<Upper>(),
                   b.block(i+s, j+block_size-1, rows-r, std::min(blocksize,cols)));
    // on right
    solve_in_place(b.block(i+s, j+block_size-1, rows-r, std::min(blocksize,cols)),
                   A.selfadjointView<Upper>());

***   The block above will work only if passed by value (temporary) but this special-case can be
handled by template specialization I suppose.

*** The free solve() functions could internally call solve_in_place() if x and b are aliases. The user will have
some choice too.


[II] x = m.lu().solve<OnTheRight | Adjoint>(b);
   seems ok to me. Is the following a little better and possible?

    x = m.adjoint().lu().solve<OnTheRight>(b);

   (see also, item [IV] below)


[III] IMHO, m.lu().solve(x)  is more readable than  m.lu().solve() * x  or m.lu().solver()*x
   - the latter is a deviation from the semantics that the inverse is an operator or function
    that operates on an input vector. If "solve" or "solver" were replaced with "inverse" this
    would be more clear. But implementation issues you refer to might make this difficult.

* I don't know how .inverse() is implemented now, but it could return an expression (proxy?) which,
  when appears in an '*' expression, could be expanded into an LU solution as opposed to multiplication
  by a computed inverse. But the Matrix class will need to have one more definition for the = operator
  to check for the inverse-expression and to evaluate the inverse from the expression template.


[IV]  Since we are applying A^{-1} or A^{-H/T}, IMHO

    b * A.adjoint().lu().inverse()

seems better to me than     b * A.lu().adjoint().inverse() (Gael)  or  b * A.lu().inverse().adjoint()  (Benoit).
This is because the first substring, A or A.adjoint(), represents the operation immediately, and then .lu().inverse()
would represent the details of the method being used to compute that inverse expression.



[V]  I think Eigen's commitment to providing a readable API is great. I think it is perfect OK to not have a "solve"
just to appease computer scientists - or to have one in addition to a more mathematically syntactic API. I am a
computer scientist myself :-)  Some good documentation should take care of concerns about in-place solution
and time-space efficiency of the algorithms.



Thanks,
Manoj