Standardize code path to reconstruct factorizations

Internals & Design Internals
,
1

What this is about

Reading the code for various factorizations in LinearAlgebra, I notized that each define a bunch of methods to reconstruct a factorization (ie obtain an ::AbstractMatrix).

Eg Cholesky has

AbstractMatrix(C::Cholesky) = C.uplo == 'U' ? C.U'C.U : C.L*C.L'
AbstractArray(C::Cholesky) = AbstractMatrix(C)
Matrix(C::Cholesky) = Array(AbstractArray(C))
Array(C::Cholesky) = Matrix(C)

while LU has

AbstractMatrix(F::LU) = (F.L * F.U)[invperm(F.p),:]
AbstractArray(F::LU) = AbstractMatrix(F)
Matrix(F::LU) = Array(AbstractArray(F))
Array(F::LU) = Matrix(F)

I am wondering if it would make sense to add default fallbacks and document which one the user would need to define for a subtype of Factorization.

Proposal

I would propose that the user defines the AbstractMatrix method, and have

AbstractArray(F::Factorization) = AbstractMatrix(F)
Matrix(F::Factorization) = Array(AbstractArray(F))
Array(F::Factorization) = Matrix(F)

This would

  1. remove repetitive code,

  2. standardize the method for future extensions, either in standard packages or user packages.

1 Like
2

Seems sensible. To be honest, I didn’t know these methods even existed.

On a tenuously related note: I’d love to have multiplication (i.e., */lmul!/rmul!) methods for <:Factorization. Multiplying with a factorized matrix is of the same cost (usually negligibly cheaper) as dividing, and this capability would mean you don’t need to keep the unfactorized version around if you require both operations. To define these methods oneself is error-prone and requires defining a new multiplication function to wrap it (or piracy).

I thought I’d mention this here since it’s very closely related to these methods. In fact, it would seem that *(::Factorization, I) (or maybe *(I, ::Factorization) would be preferable, depending on which side any pivoting occurs) would generalize the above AbstractMatrix methods (although is not necessarily how they should be implemented).

It would seem that something more along the lines of

Matrix(F::Factorization) = convert(Matrix, AbstractArray(F))

would save a copy in many cases, since the AbstractArray method would sometimes return a Matrix already.