Extending Linear Algebra for Units

Internals & Design
1

I’m the author of the FlexUnits.jl library and I’ve currently built some functionality around linear algebra and unitful matrices. The key behind efficient linear algebra is to factor the units of an MxN matrix into a LinmapQuant containing a pure numerical matrix and a factored-out unit matrix.

LinmapQuant(m::AbstractMatrix{Real}, dims::DimsMap)

Where DimsMap contains SI units, with a canonical representation as:

DimsMap(u_fac, u_in, u_out)

where you can reconstruct a matrix of units with

U = u_fac * u_out * inv.(u_in)'

This has N+M-1 degrees of freedom (the first element of u_in and u_out is always dimensionless). This allows resolving units to be operations of order O(N+M), and in some cases O(N) or even O(1).

My question is optimal design around factorizations. My current experimental design is to build a wrapper around a factorization using FactorQuant{F} <: Factorization

struct FactorQuant{F, D<:AbstractDimLike, U<:AbstractDimsMap{D}}
    factor :: F
    dims  :: U 
end

which does not subtype to anything. It separates the units from the factorization and mirrors the API of an internal factorization. The advantage of this approach is that

  1. It avoids ambiguity
  2. It factors out units in a way that is usually more efficient than the underlying factorization

The problem is that methods like cholesky(m::LinmapQuant) will not produce a Cholesky decomposition, but instead a FactorQuant{Cholesky} which behaves like a Cholesky decomposition. You can call ch.L and ch.U on it because it forwards those methods to the underlying numerical Cholesky and wraps the appropriate units around the response.

Would it be more proper to actually produce a “Cholesky{LinmapQuant}” and overload the “getproperty” function to produce LinmapQuant instead? From what I can see the pros and cons of this general approach would be:

Pros:

  • Type-Consistency: cholesky(m::LinmapQuant) returns an actual Cholesky factorization, not FactorQuant{Cholesky}
  • Integration: Many algorithms subtype onto a Cholesky, or Eigen decomposition

Cons:

  • Performance (potentially small): Moving a LinmapQuant inside a factorization hides information about the unit structure. It’s always better to resolve the units of a matrix/factorization separately due to degrees of freedom. This may be mitigated by specializing on f(m::LinmapQuant) for constructors and getindex(fac::F{LinmapQuant})
  • Ambiguity (potentially massive): Specializing on LinmapQuant can lead to ambiguities with multiple arguments. Pointing to specializations on inner arguments is tedious and manual, the FactorQuant does this automatically.

This also gets into the issue of, what should “LowerTriangular(m::LinmapQuant)” return? LowerTriangular{LinmapQuant} is consistent, but LinmapQuant{LowerTriangular} is more efficient. This is essentially the problem of a factorization of a factorization, where one factorization is more efficient. I’m wondering if the community knows which strategy is superior and why.

2 Likes