Re: [eigen] solve API |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] solve API*From*: Manoj Rajagopalan <rmanoj@xxxxxxxxx>*Date*: Mon, 28 Jun 2010 10:42:19 -0400*Organization*: EECS Dept., University of Michigan, Ann Arbor, MI, USA

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

**Follow-Ups**:**Re: [eigen] solve API***From:*FMDSPAM

**Re: [eigen] solve API***From:*Gael Guennebaud

**References**:**[eigen] solve API***From:*Gael Guennebaud

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] sparse colwise** - Next by Date:
**[eigen] documentation: linear algebra/decompositions** - Previous by thread:
**Re: [eigen] solve API** - Next by thread:
**Re: [eigen] solve API**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |