Eye in julia 0.7

I think I found a replacement for eye in julia 7.0

julia> c=rand(Int8,4,4)
4×4 Array{Int8,2}:
   44   -79   36  116
  -14   -64  -49  -52
   77  -110  -20  118
 -119   -57  -78   90

julia> @btime Matrix{eltype($c)}(I,size($c))
  47.614 ns (1 allocation: 112 bytes)
4×4 Array{Int8,2}:
 1  0  0  0
 0  1  0  0
 0  0  1  0
 0  0  0  1

julia> @btime $c^0
  33.817 ns (1 allocation: 112 bytes)
4×4 Array{Int8,2}:
 1  0  0  0
 0  1  0  0
 0  0  1  0
 0  0  0  1

Am I right? Is c^0 a good replacement for eye(c) for square matrices?

I do find it annoying that the constructor for a properly typed I in 0.7 is so verbose.

I think it is intended that you try to use I itself wherever possible, rather than constructing a new matrix, and the thinking was probably that the latter case is rare enough that it doesn’t need condensed syntax.

Anyway, I’m seeing that for larger matrices c^0 is actually (very slightly) slower than the Matrix constructor. In cases like this where the performance differences are very slight it depends to be best to use the intended functionality, though I’m afraid I can’t comment further on exactly what c^0 is doing.

1 Like

c^0 just boils down to a call to one(c). It’s slightly different than Matrix{eltype(c)}(I, size(c)) in that it returns the multiplicative identity — which might be different for unitful quantities.

If one(x) is faster than Matrix{T}(I, size(x)), then that’s a performance bug in the Matrix constructor.

1 Like

You can just use one(A) or oneunit(A) (depending on how you want to handle the dimensionful case) for the equivalent of eye(A) for square matrices.

(eye also worked for non-square matrices, but I don’t know of a non-verbose equivalent on 0.7.)

4 Likes

with the same example as above

julia> @btime one($c)
  27.680 ns (1 allocation: 112 bytes)

quite a bit faster than the Matrix constructor…

The cost of either is trivial in itself, and is dwarfed by anything meaningful you may do with the matrix, eg

using BenchmarkTools
c = rand(Int8, 10, 10)
altone(A) = Matrix{eltype(A)}(I,size(A))
twice(A) = one(A) * 2
alttwice(A) = altone(A) * 2

gives

julia> @btime twice($c);
  142.838 ns (2 allocations: 1.06 KiB)

julia> @btime alttwice($c);
  157.760 ns (2 allocations: 1.06 KiB)

Agreed, but the question is: when I have a square matrix A, should I type the suggested replacement Matrix{eltype(A)}(I,size(A)) for eye(A) while one(A) seems to give the same result and is faster?

My understanding was that the suggested solution is one, so it is not clear to me why one would use something else.

julia> eye(rand(4,4))
WARNING: Base.eye is deprecated: it has been moved to the standard library package `LinearAlgebra`.
Add `using LinearAlgebra` to your imports.
 in module Main
┌ Warning: `eye(A::AbstractMatrix{T}) where T` has been deprecated in favor of `I` and `Matrix` constructors. For a direct replacement, consider `Matrix{eltype(A)}(I, size(A))`.If `eltype(A)` element type is not necessary, consider the shorter `Matrix(I, size(A))` (with default `eltype(I)` `Bool`).
│   caller = top-level scope
└ @ Core :0
4×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0
2 Likes

I agree that suggesting one(A) would be a better choice there.

1 Like

No, because one only works for square matrices.

Just out of curiosity, what’s the motivation for deprecating eye?

Certainly, but asking for a multiplicative identity is a probably the most common case.

See make I(n) behave as `Diagonal` / `eye`? · Issue #23156 · JuliaLang/julia · GitHub and the linked issues, PR’s

2 Likes

I’ve always wondered why there isn’t something like

struct Identity{T} <: AbstractMatrix{T}
    n::Int
    m::Int # perhaps
end

i.e., a dedicated isbits identity matrix type that is an AbstractMatrix, and has both size and eltype. This type would also be useful for cases such as https://github.com/JuliaLang/julia/issues/23156#issuecomment-340216693.

I understand that part of the objective of the change was to reduce the number of ways to create an identity matrix, but I find the typelessness and shapelessness of I only useful in a subset of cases, and the Diagonal version allocates (unless you make a special isbits AbstractVector type for the diagonal).

1 Like

I’d be curious to hear what cases you have that require shape but not mutability.

The deprecation has to give code that always works, not just code that works in the common case.

Well, https://github.com/JuliaLang/julia/issues/23156#issuecomment-340216693 for one. It could also be used as an implicit shape check in stuff like

I33 = Identity{Float64}(3, 3)
[rand(4, 4) I33;
 I33 rand(4, 4)]

This would work if I33 were replaced with I, but there may be cases where that behavior is actually not desired.

1 Like

Should not the deprecation help the programmer as much as possible? Here the most common subcase
has a nicer/shorter deprecation, so should not this be mentioned in the deprecation message?

1 Like

Auto-deprecations are done in Julia with @deprecate oldfunc(args...) newfunc(newargs...), where newfunc should give the same result as oldfunc in all cases.

It is possible to make custom deprecation messages, but it is more work and we usually don’t bother (since the deprecation message will go away entirely in the subsequent release anyway).

Note also that, in this particular case, you are really strongly encouraged to re-think constructing a dense identity matrix entirely — there is almost always a better way to do things. What did you need eye for anyway?

1 Like