I need an isinvertible
function, where isinvertible(x)
is true iff inv(x)
is defined, and isone(inv(x) * x)
. That’s to rule out 0
, since inv(0)
is defined, but returns Inf
.
In the Julia LU decomposition, There is a test for zero. That is to protect the taking of inv
a few lines later. So that iszero
test really wants to be testing to see if the element is invertible. This matters for me, because I want to take inverses of matrices of modular integers, and often there are non-zero elements which are not invertible.
isinvertible(x) = applicable(inv, x) && isone(inv(x) * x)
1 Like
The condition that the result be one is not necessarily true, eg in interval arithmetic.
1 Like
@bobcassels has another question outstanding: How to test whether comparison is defined for a type?
I thought that the `applicable function might help there too, but it doesn’t.
julia> applicable(<, Mod29(1) , Mod29(3))
true
but this fails
julia> Mod29(1) < Mod29(3)
ERROR: < not defined for Mod29
Because
(< )(x::T, y::T) where {T<:Real} = no_op_err("<" , T)
# see base/promotion.jl line 426
That is, the <
method applies to the arguments, but unfortunately that method forces an error.
So my suggestion to use the applicable
function is not foolproof
What kind of object you are expecting as input? Number? Matrices?
Sorry, just saw the end of the post, it was cropped on my tablet.
The problem you saw can be avoided by a try ... catch
block.
In my case, I am interested in anything that might be an element in a matrix passed to LU decomposition. For now, my modular integers, <:Real.
Of course I can use try ... catch
. But that doesn’t seem “in the spirit of Julia”, at least as I understand that spirit. I want to put this in the middle of the standard library LU decomposition function.
using LinearAlgebra
isinvertible(x) = issuccess(lu(x, check=false))
2 Likes
julia> isinvertible(3)
ERROR: UndefVarError: lu not defined
using LinearAlgebra
first. I will edit the comment.
I’m not trying to check the result of lu. I want to change the source of lu to check to see if an element is invertible, and do something different if it is not.
Can you give examples with inputs and desired outputs of this isinvertible
function?
Imagine I have defined a type Mod26 that is integers mod 26 [where Mod26(x) is invertible iff gcd(x, 26) == 1]
isinvertible(100) # true
isinvertible(0) # false
isinvertible(-1.5) # true
isinvertible(-0.0) # false
isinvertible(Mod26(2)) # false
isinvertible(Mod26(7)) # true
isinvertible(Rational(1) + Rational(3)im) # true
isinvertible(0+0im) # false
These should work:
isinvertible(x) = !iszero(x)
isinvertible(x::Mod26) = gcd(x, 26) == 1
Let’s pause this discussion for now. It’s not going in a useful direction.
Yes, I can define it myself. But I want to change the system’s LU decomposition, so that it works for more cases. So I want this isinvertible
function to be part of the core, so it can be used by the standard library. I will make a pull request for the changes I want to make to Julia, and then we can discuss this in that context.
I want the standard library’s matrix inversion to work on matrices of modular numbers. It does not currently, for the reasons I’ve pointed out. I am proposing changes to allow it to work. As I said, I will propose that change, and we can discuss in that context.
If I understand your situation correctly, it sounds like this should be solved by multiple dispatch + extension:
- Add an
isinvertible
function to Base that covers all relevant default types
- Change the LU source from
!iszero
to isinvertible
- Extend
isinvertible
in your code to work with your type:
import Base.isinvertible
isinvertible(x::Mod26) = gcd(x, 26) == 1
Btw, if you don’t want to go through the process of changing Julia itself, perhaps you could just extend iszero(::Mod26)
? A bit of a hack, but it might work for you.
2 Likes
Please don’t suggest this Exceptions should be reserved for exceptional circumstances, not situations like this.
3 Likes
Note that determining if an object in a nontrivial algebra is (numerically) invertible often has a cost that is comparable to inverting it. Checking if something is invertible without reusing the intermediate representation (eg LU for matrices, etc) may be inefficient.
It may be better to just explain your problem in detail, maybe it has a better solution. Eg if you want to generalize LU, then do that instead.
1 Like