# Extending Linear Algebra for Units

**URL:** <https://discourse.julialang.org/t/extending-linear-algebra-for-units/139761>\
**Category:** Internals & Design\
**Created:** [September 30, 2026, 3:53pm UTC](https://discourse.julialang.org/t/extending-linear-algebra-for-units/139761 "2026-09-30T15:53:29Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Deduction42](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/deduction42/32/9206_2.png) [@Deduction42](https://discourse.julialang.org/u/Deduction42)\
**Post date:** [September 30, 2026, 3:53pm UTC](https://discourse.julialang.org/t/extending-linear-algebra-for-units/139761/1 "2026-09-30T15:53:29Z")

</div>

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.

```julia-auto
LinmapQuant(m::AbstractMatrix{Real}, dims::DimsMap)

```

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

```julia-auto
DimsMap(u_fac, u_in, u_out)

```

where you can reconstruct a matrix of units with

```julia-auto
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`

```julia-auto
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.