Similar identity matrix

General Usage
1

What is the idiomatic way to construct a full m×m or n×n identity matrix similar to a given m×n matrix A (static or not)?

one(A*A') or one(A'*A) does the trick but is of course not what I want.

2

Do you really need one? Can’t you just sidestep the question with LinearAlgebra.I?

3

Yes, it is the initial value of a matrix-valued ODE and should be of the same type as later values.

4

One way is

julia> using LinearAlgebra

julia> A = rand(4,5);

julia> diagm(0=>fill(1., size(A,1)))
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

julia> diagm(0=>fill(1., size(A,2)))
5×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0

or even diagm(0=>fill(one(eltype(A)), size(A,1))).

Note that in Julia 0.6 you could do

julia> eye(size(A,1))
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

julia> eye(size(A,2))
5×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0
5

Another way would be Matrix(one(eltype(A))I, size(A,1), size(A,1)) I guess.

6

Thank you, that’s it (almost)!

julia> A
5×4 Array{Float64,2}:
 0.896177   0.0597406  0.611605   0.231433
 0.119473   0.998269   0.0299731  0.947813
 0.0465002  0.431471   0.321885   0.20422 
 0.641878   0.0221313  0.0604966  0.517904
 0.958566   0.675325   0.0570553  0.603291

julia> typeof(A)(one(eltype(A))I, size(A,1), size(A,1))
5×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0
7

Oh, I see.

8

For static arrays it is unfortunately slightly different

julia> similar_type(B, Size(size(B,1),size(B,1)))(I)
5×5 SArray{Tuple{5,5},Float64,2,25}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0
9 
julia> using LinearAlgebra

julia> A = rand(5,4);

julia> copyto!(similar(A, size(A,1), size(A,1)), I)
5×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  1.0
1 Like
10

I don’t think this works for static matrixes:

julia> using LinearAlgebra, StaticArrays
julia> A = @SMatrix rand(5,4);
julia> typeof(copyto!(similar(A, size(A,1), size(A,1)), I))
Array{Float64,2}

because the generic copyto! is used.

This is could be a bug in StaticArrays, but at the same time I suspect that the design of the library which needs this could be improved, well-designed Julia code should either not require tricks like this, or provide helper functions for it.

11

For this use-case it would be nice provide a second argument like Val(:Left) to one which allows to specify to take the left or a right multiplicative unit if they are not equal.

1 Like