RFC: Improved tridiagonal and new Woodbury matrix algebra

[timholy]

This commit:

• Converts Tridiagonal to being Lapack-compatible
• Wraps Lapack's major tridiagonal routines
• Adds a new "Woodbury" type for solving equations using the Woodbury matrix identity
• Provides a set of tests for this functionality

Issues that might be worth discussing before merging this:

• Should the tridiagonal routines be moved into linalg_lapack.jl and friends? Conceptually, tridiagonal routines are sparse, which is why I've been putting them in extras/sparse.jl. However,
  • Inside linalg_lapack.jl, at least some "banded" matrix routines, like gbtrf, are Julia-wrapped. Tridiagonal matrices are simply an extreme case of banded matrices. (As far as I can tell, there's no way to actually use these banded matrix routines currently, because the basic Lapack wrappers are not exported and there is no higher-level routine that calls them.)
  • There's some awkwardness to re-loading the OpenBLAS library, particular if the user has chosen a different BLAS library
• Currently I don't provide an lu wrapper for LU-decomposition, but there is a "private" (underscore) wrapper for tridiagonal LU decomposition, a wrapper to use the LU decomposition to solve equations, and a test of this functionality in test/sparse.jl. The reason I don't provide a nice lu function is that the LU-decomposition of a tridiagonal matrix requires the second off-diagonal. This can actually be stored in the Tridiagonal type (and I do so), but the meaning of the object is not what it once was. There are several possible ways around this:
  • Have lu(M::Tridiagonal) return full (i.e., dense) L and U matrices. This stinks so bad, it is hardly worth even mentioning.
  • Have lu(M::Tridiagonal) do what the underscore Lapack wrapper does now: return an object of type Tridiagonal, but this is something of a lie, because the second off-diagonal stored in the dutmp field is actually meaningful (it just happens to be convenient storage that is already allocated). The second output would be ipiv or a permutation constructed from it. In other words, the return-value syntax of lu would differ from that of dense LU-decomposition.
  • Create BandedMatrix types in Julia. This may be the most ambitious route, with the advantage that the current return syntax of lu could be preserved. However, there is a performance penalty here, because the LU decomposition of a tridiagonal can be stored and manipulated more efficiently than what is required for general banded matrix support.
  • Introduce a new alternative syntax lu(LU::LUMatrixDecomposition, A::Matrix), a function solve(LU::LUMatrixDecomposition, rhs::StridedVecOrMat), and types to support this. This syntax allows you to control the output format of the result, and uses pre-allocated storage consistent with issue Functions that produce results in pre-defined outputs #1115. It would also allow for maximally-efficient usage of the Lapack routines, without the need to call tril and triu on the result. That said, compared to the O(N^3) cost of LU-decomposition, the O(N^2) overhead of tril and triu should be pretty inconsequential. So the main advantage of this route is simply that you can customize the return type appropriately to the input type, which is something you can't do with the current lu syntax.

For this second issue, I personally think the fourth option is the best. Interested in hearing feedback from others.

[dmbates]

Thank you for converting the Tridiagonal type to be Lapack-compatible and for using the Lapack tridiagonal routines. I think it is definitely a good practice to use Lapack capabilities, when available, not only for consistency of approach but also because there may be behind-the-scenes work to optimize them.

I would encourage putting the Lapack wrappers in linalg_lapack.jl I think the plan is to have that file be the place where all the Lapack wrappers are defined and put them inside a module. Most user code won't access those wrappers directly but other modules can import the namespace and access the wrappers through a qualified name. There really should only be one _jl_liblapack and one _jl_libblas handle on Lapack and/or BLAS libraries, and those are the ones created in start_image.jl

I'm not sure I would put the tridiagonal code in *sparse.jl files. To me the *sparse.jl files contain code that handles triplet and CSC (compressed sparse column) representations of matrices. Even though banded matrices are sparse they are patterned and do not require the (i, j, x)-like representation of the sparsity. However, that is just my opinion.

Have you considered creating a factorization type rather than trying to create methods for the lu() function? Your fourth point, "Introduce a new alternative syntax" is essentially what is implemented in base/factorizations.jl, isn't it? Unfortunately factorizatilons.jl defines an LU type as the LU decomposition of a dense general matrix. I suppose it would be possible to create an abstract LU type with concrete subtypes for LU of a general matrix and LU of a tridiagonal matrix or just adopt the convention the a factorization type like LU or Cholesky or QR is the decomposition of that type for a general matrix and decompositions for patterned matrix types would have other names.

[timholy]

I had missed base/factorizations.jl! Yes, that's essentially exactly what I'm proposing. As you say, we'd need a slightly more flexible implementation, but I think that's merely details at this point.

One point of confusion (and why I missed it) is because lu is defined in base/linalg_lapack.jl. Coming from a Matlab heritage, I found lu and looked no further. (I also note that it's the only variant that is mentioned in the documentation.) Not entirely sure what to do about this. Perhaps the best bet would be to make lu produce an error when called with Tridiagonal inputs, and force people to use the LU syntax. The error message might even be created to be helpful, and direct people appropriately :-).

I think you've answered both of my major issues: move Tridiagonal into base, and make LU-decomposition work like base/factorizations.jl. I'll give a little more time for others to chime in, too, in case there's any disagreement, but I agree with both recommendations.

This does leave a smaller issue: given how you're thinking about extras/sparse.jl, Woodbury also doesn't fit. Should we put it wherever Tridiagonal ends up? It doesn't use Lapack, of course.

[timholy]

@dmbates , see if you like this better.

If it is acceptable, I guess the only remaining question is whether to merge as-is, or if I should squash it and merge directly into master (to reduce repository churn).

[dmbates]

After these changes it looks fine to me. You may want to squash the multiple changes or, if you are as inept with git as I am, it may not be worth the effort.

It might be good to have @ViralBShah or @StefanKarpinski comment before it gets committed.

[ViralBShah]

This looks good. Would be nice to squash, since we don't really need the history for the modifications to sparse.jl here. @timholy, can you rebase and merge?

[ViralBShah]

Also, we do need to get rid of lu returning matlab style decomposition matrices, and make it simply give an LU matrix factorization type. We need to use the factorization style interface throughout all our dense linear algebra.

After that, I think we need to move all the blas calls into a BLAS module, lapack calls into an LAPACK module, and all the linear algebra stuff into a Linalg module. Then, all of linalg should be loaded by default on startup, but the rest of the stuff can stay hidden. I guess we need to improve docs as well.

Not saying all of this should be done in this pull request, but as the way forward so that our linear algebra stuff is release-ready.

[timholy]

Yes, I can squash. Figured out how to do that a couple of days ago. Slowly getting more sophisticated at git...

As you say, I will not tackle the module stuff as part of this. I really do like the idea of getting rid of the matlab style decomposition, I will do that gladly. I'd argue we should still have an lu (i.e., lowercase) function, which can return LU types.

[ViralBShah]

Yes, we should retain the matlab names, but have them return Factorization objects for lu, qr, chol, svd, etc. These Factorization objects can all then have a common set of functions like solve, and some extensions for specific cases.

[toivoh]

Extensions including extracting e g the actual L and U factors? (and permutation P).
Sometimes you do need the factors themselves.

[ViralBShah]

I was thinking about an extract or factors method for factorization objects? So, if I want the factors, I do:

(l,u,p) = factors(lu(A))

-viral

On 03-Sep-2012, at 12:46 PM, toivoh wrote:

Extensions including extracting e g the actual L and U factors? (and permutation P).
Sometimes you do need the factors themselves.

—
Reply to this email directly or view it on GitHub.

[toivoh]

Looks good to me.

[dmbates]

I was going to suggest using extract' orfactors' myself. In the Matrix package for R we used extract' butfactors' applied to a Factorization object probably makes more sense.

There may be some confusion if the Matlab-compatible names like chol and lu do not return Matlab-compatible values. Why not make lu(A) return factors(LU(A)) and in the documentation state that direct use of lu() is deprecated?

[ViralBShah]

I was actually suggesting that lu(A) return a type LUFactor (I presumed LU was a type that was short for LUFactor). I am not sure of having both lu and LU, both having different meanings.

Do you think it would be a major issue with matlab? After all, they are returning factorizations, which is what matlab does as well. I feel it is reasonable to retain matlab compatible names, but return factorization objects. After all, our factorization objects do support the use of \ which is what one would do in matlab too.

-viral

On 03-Sep-2012, at 8:43 PM, dmbates wrote:

I was going to suggest using extract' orfactors' myself. In the Matrix package for R we used extract' butfactors' applied to a Factorization object probably makes more sense.

There may be some confusion if the Matlab-compatible names like chol and lu do not return Matlab-compatible values. Why not make lu(A) return factors(LU(A)) and in the documentation state that direct use of lu() is deprecated?

—
Reply to this email directly or view it on GitHub.

[StefanKarpinski]

For what it's worth, I like calling the type LU, especially if all the other ones follow that pattern.

[ViralBShah]

It is not nice to keep using up short names in libraries. LU is something that I would love to use in my user code.

-viral

On 04-Sep-2012, at 2:10 AM, Stefan Karpinski notifications@github.com wrote:

For what it's worth, I like calling the type LU, especially if all the other ones follow that pattern.

—
Reply to this email directly or view it on GitHub.

[StefanKarpinski]

I guess, but we don't want to become like Mathematica — I do not want to type SingularValueDecomposition[X].

[timholy]

Largely consistent with the advice here so far, I propose that

• Capitals be used to refer to types, e.g., LUTridiagonal, LUBand (not yet implemented), and either LU or LUDense. I don't have strong feelings about the latter, happy to accept guidance.
• lu be the generic function that would compute the appropriate LU-decomposition for any matrix type. Its return type would be predictable from its inputs but would be of varying types depending on the type of the input, i.e., LUTridiagonal if the input is Tridiagonal.

The second point explicitly breaks Matlab compatibility. I am personally in favor of this particular breakage, and like the idea of factor as a separate function. On balance I think this is a pretty reasonable breakage because it results in an error rather than a bug. The biggest negative is the nature of the error message:

julia> L, U, P = LU(m)
type error: tupleref: expected (Any...,), got LU{Float64}

(I used current naming just to provide a real-world example, it would change to lu according to this proposal.) I suspect that a newcomer to Julia would find this error much more confusing than something like "lu() returns a LU-factorization type, use only one output. To extract the individual matrices, use factor() on the output." But without an nargout, I don't think there's a way to provide this type of error message.

[ViralBShah]

LUDense is the right compromise for now. Rest sounds great.

-viral

On 04-Sep-2012, at 6:04 AM, Tim Holy notifications@github.com wrote:

Largely consistent with the advice here so far, I propose that

Capitals be used to refer to types, e.g., LUTridiagonal, LUBand (not yet implemented), and either LU or LUDense. I don't have strong feelings about the latter, happy to accept guidance.
lu be the generic function that would compute the appropriate LU-decomposition for any matrix type. Its return type (see next bullet) would be predictable from its inputs but would be of varying types depending on the type of the input, i.e., LUTridiagonal if the input is Tridiagonal.
The second point explicitly breaks Matlab compatibility. I am personally in favor of this particular breakage, and like the idea of factor as a separate function. On balance I think this is a pretty reasonable breakage because it results in an error rather than a bug. The biggest negative is the nature of the error message:

julia> L, U, P = LU(m)
type error: tupleref: expected (Any...,), got LU{Float64}
(I used current naming just to provide a real-world example, it would change to lu according to this proposal.) I suspect that a newcomer to Julia would find this error much more confusing than something like "lu() returns a LU-factorization type, use only one output. To extract the individual matrices, use factor() on the output." But without an nargout, I don't think there's a way to provide this type of error message.

—
Reply to this email directly or view it on GitHub.

[ViralBShah]

Mathematica is overboard with names. Square brackets for function calls is also a bit weird for me, but I am digressing.

-viral

On 04-Sep-2012, at 5:50 AM, Stefan Karpinski notifications@github.com wrote:

I guess, but we don't want to become like Mathematica — I do not want to type SingularValueDecomposition[X].

—
Reply to this email directly or view it on GitHub.

[timholy]

Closed with commits 826c42f and cbaad36 to master.