When to throw errors when calling inv(::AbstractCliffordNumber)

[brainandforce]

The function inv(::AbstractCliffordNumber) is implemented using the internal function CliffordNumber.versor_inverse, which just calculates x' / x^2, with no regard to whether the result actually approximates the inverse reasonably well. The overload of Base.inv adds logic to check if the result is actually an inverse, and throws a CliffordNumbers.InverseException for Clifford numbers without an inverse.

One example of a CliffordNumber without an inverse is CliffordNumber{VGA(2)}(1, 1, 0, 0), and this inverse call throws an error:

julia> x = CliffordNumber{VGA(2)}(1, 1, 0, 0)
4-element CliffordNumber{VGA(2), Int64}:
1 + 1e₁

julia> inv(x)
ERROR: CliffordNumbers.InverseException: The result of x' / abs2(x) was not an inverse.
Note that not all Clifford numbers will have an inverse.
Stacktrace:
 [1] inv(x::CliffordNumber{VGA(2), Int64, 4})
   @ CliffordNumbers ~/git/CliffordNumbers.jl/src/math/inverse.jl:43
 [2] top-level scope
   @ REPL[159]:1

julia> x * CliffordNumbers.versor_inverse(x)
8-element CliffordNumber{VGA(3), Float64}:
1.5 + 0.5e₁

In positive-definite metrics, all 0-blades, 1-blades, and their complements will have inverses, but this is not the case in other metrics. I wrote an optimization to handle these cases, but this causes null vectors in other metrics to return an inverse consisting of Inf or NaN entries:

julia> inv(KVector{1,STA}(1, 1, 0, 0))
4-element KVector{1, LGAWest(3), Float64}:
Infγ₀ + Infγ₁ - NaNγ₂ - NaNγ₃

The question here is: should this behavior be the default, or should we through a CliffordNumbers.InverseException for these cases?

[moble]

A couple relevant data points:

julia> inv(0.0)
Inf

julia> inv(0.0 + 0.0im)
NaN + NaN*im

FWIW, my opinion is that anyone sophisticated enough to be using these GAs should be sophisticated enough to understand these cases, and must be willing to avoid or handle them. I wouldn't expect inv to ever throw.

[brainandforce]

I wouldn't expect inv to ever throw.

Would you also expect this to be the case for idempotents/projectors/other elements that cannot have an inverse? The reason I ask is because if you try to invert a singular matrix, it throws an error (version 1.10.5), fresh REPL:

julia> inv([1 0; 0 0])
ERROR: LinearAlgebra.LAPACKException(2)
Stacktrace:
 [1] chklapackerror
   @ /usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:40 [inlined]
 [2] trtrs!(uplo::Char, trans::Char, diag::Char, A::Matrix{Float64}, B::Matrix{Float64})
   @ LinearAlgebra.LAPACK /usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:3564
 [3] generic_trimatdiv!
   @ /usr/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:830 [inlined]
 [4] _ldiv!
   @ /usr/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:752 [inlined]
 [5] ldiv!
   @ /usr/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:757 [inlined]
 [6] inv(A::LinearAlgebra.UpperTriangular{Int64, Matrix{Int64}})
   @ LinearAlgebra /usr/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:797
 [7] inv(A::Matrix{Int64})
   @ LinearAlgebra /usr/share/julia/stdlib/v1.10/LinearAlgebra/src/dense.jl:913
 [8] top-level scope
   @ REPL[4]:1

The semantics of the package intentionally treat multivectors as a numeric type (subtyping Number), but there are some areas where the properties diverge (for instance, OddCliffordNumber(::CliffordNumber{Q,T}) performs a grade projection, but convert(::Type{<:OddCliffordNumber}, ::CliffordNumber{Q,T}) performs conversion and throws an InexactError if OddCliffordNumber would strip even grade coefficients (for most Number subtypes, these behave identically).

I'm personally okay with inv throwing an exception for some cases, but I also want consider your perspective here too and minimize the need for exception handling. Perhaps it would be useful to supply something along the lines of checked_inv that can perform this check if explicitly specified, and requiring that inv never throws, but this would be a breaking change for 0.2.0.

There may also be some notion of pseudoinverse that would work here - perhaps the simple pinv(x) = x' / abs2(x), perhaps something more robust. I also know there are algorithms for performing inversions of arbitrary multivectors in algebras of low dimension, but that isn't something I've taken the time to implement yet - that would be a separate issue/PR.

[moble]

if you try to invert a singular matrix, it throws an error

That's a good point. I guess I should have clarified that I don't have a strong opinion on this. And even by my own logic, the sophisticated user should be just as capable of avoiding and handling exceptions as infs and nans.

In this sort of situation, I think it pays to back up and think about fundamental priorities. I suppose correctness should be first, but you've probably got that in the bag. So maybe the second should be speed. I suppose checking for error conditions would be slow for built-in floats and complex, but as you get to more complicated structures (like matrices and multivectors) maybe the marginal cost of a check isn't much, so you get to favoring errors.

Maybe the next priority should be simplicity of your code. Depending on exactly what's happening, I suppose that might favor nan/inf. The only relevant issue here that I've run into is ensuring that the code's structure is simple enough to be differentiated correctly near these special values.

Otherwise, you're the author, so if there's no great argument one way or another, just do what you feel is right, but be sure to document it! :)