search

gonum / matrix

Matrix packages for the Go language [DEPRECATED]
446 stars 51 forks source link

Redesign Solve #169

Closed btracey closed 9 years ago

btracey commented 9 years ago

(Adapted from an off-line discussion)

Lapack is almost at a state where it can implement the full functionality of the current Solve.

I don’t think the current Solve function is idiomatic within mat64. The result of the operation is a matrix, so it should be stored in-place into a type like any other operation. I have found that linear solving for a vector is very common, so I would suggest

func (m *Dense) Solve(a, b Matrix) error // possibly bool
func (v *Vector) SolveVec(a Matrix, b *Vector) error // possibly bool

I have learned in dealing with Lapack that the factorizations have special forms for storing data. The PLU decomposition stores L and U in a single matrix, as does the QR factorization. I believe this merits the continued existence of the LUFactors and QRFactor types. I would guess we can find better names (at the very least we should pick one of Factor and Factors), but the current ones aren’t terrible. Possibly we could have the name just be “LU”, and then we could have 

// Factorize finds the LU factorization of the matrix
func (lu *LU) Factorize(a Matrix)

As now, the user can extract the types. It possibly should be something like

// GetL extracts the L of an LU decomposition
func (t *Triangular) GetL(lu *LU)

It is common to want to perform Solve many times using the same matrix. It is much more efficient to perform the factorization once and then use the factorization many times to perform different solves. Additionally, the various decompositions have differing tradeoffs when it comes to linear solving. For example, LU has a lower constant, but is less robust. It seems nice to enable solving using specific methods for users with special needs (or, if Solve fails, they can error check and try QR or whatever). Both of these imply having a special solve like we now have with Cholesky

// SolveLU computes the solution to a set of equations using a pre-computed LU factorization. 
func (m *Dense) SolveLU(lu *LU, trans bool, b Matrix)

The trans field is necessary because LU doesn’t actually implement Matrix. We would similarly provide a method for QR and the corresponding *Vector.

Lastly, I feel that the LU/QR decompositions merit certain methods. I have been working on an implementation of the Simplex algorithm for solving linear programs. The core of the simplex algorithm is a set of linear solves. First, there are 3 or 4 linear solves per iteration that all use the same matrix inverse (hence the first suggestion). Second, each iteration of the simplex algorithm performs a rank-one update to the inverted matrix. It is evidently possibly to apply the rank one update to the factorization directly. Doing this drops the inner loop of the simplex from m^3 to m^2, which I feel is significant enough to add a RankOne() method on to the factorization. Methods would of course be added only as they are needed.

kortschak commented 9 years ago

SGTM to reflect result of off-line discussion.