Open brainandforce opened 6 days ago
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.
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.
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! :)
The function
inv(::AbstractCliffordNumber)
is implemented using the internal functionCliffordNumber.versor_inverse
, which just calculatesx' / x^2
, with no regard to whether the result actually approximates the inverse reasonably well. The overload ofBase.inv
adds logic to check if the result is actually an inverse, and throws aCliffordNumbers.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: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
orNaN
entries:The question here is: should this behavior be the default, or should we through a
CliffordNumbers.InverseException
for these cases?