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

It’s not clear how much you will gain from this, since if a function was specifically written for a Cholesky decomposition, there’s a good chance it’s not expecting the factors to be dimensionful.

3

if a function was specifically written for a Cholesky decomposition, there’s a good chance it’s not expecting the factors to be dimensionful.

I’m no sure about that if the Cholesky decomposition is performed on a LinmapQuant which is an AbstractArray. I work in statistics and simulation of dynamic systems, and cholesky factorizations often deal with coveriance matrices which have symmetric units. I would expect a cholesky decomposition on such an object to result in dimensionful factors. You could argue that eigendecomposition for PCA should be dimensionless because it’s ideally performed on correlation matrices. However, analyzing dynamic system stability would involve taking eigenvalues of a transition matrix which has (idempotent) units.

This is a hard choice. The FactorQuant approach is definitely the easier, safer option and results in fewer invalidations. It also a clean level of abstraction that makes sure that the unit structure factorization is prioritized. The parametric approach (Cholesky{LinmapQuant}) leverages more of Julia’s potential, enabling you to insert dimensionful operations into other functions out of the box, but the maintainers won’t be expecting your types and won’t be testing against them. Something that works one day could be broken the next unless you can convince the maintainers to make it work for FlexUnits.jl which is what I think @stevengj is getting at. Parametric multiple dispatch is insanely powerful but very brittle if people aren’t pulling together.