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.
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.
Do you really need one? Can’t you just sidestep the question with LinearAlgebra.I
?
Yes, it is the initial value of a matrix-valued ODE and should be of the same type as later values.
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
Another way would be Matrix(one(eltype(A))I, size(A,1), size(A,1))
I guess.
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
Oh, I see.
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
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
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.
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.