You can use Diagonal?
3 Likes
OK – for a vector v, diagm(v) and Diagonal(v) give the same result? Or does Diagonal(v) avoid storing off-diagonal elements?
Wow!
julia> M1 = Matrix{Integer}(I,20,20);
julia> M2 = diagm(fill(1,20));
julia> M3 = Diagonal(fill(1,20));
julia> sizeof(M1), sizeof(M2), sizeof(M3)
(3200, 3200, 8)
It is possible to use
A = I(2)
2×2 Diagonal{Bool,Array{Bool,1}}:
1 ⋅
⋅ 1
4 Likes
Superb! So both of my questions are answered, I guess:
julia> using LinearAlgebra, Random;
julia> n = 4;
julia> cp = randperm(n);
julia> A = rand(-9:9,n,n)
4×4 Array{Int64,2}:
7 6 -9 5
-1 5 7 -8
-3 -2 0 -3
3 4 -4 5
julia> P = I(n)[:,cp];
julia> all(A[:,cp] - A*P .== 0)
true
julia> rp = P*(1:n);
julia> all(A[rp,:] - P*A .== 0)
true
There is an I(n) constructor, eg
julia> I(3)
3×3 Diagonal{Bool,Array{Bool,1}}:
1 ⋅ ⋅
⋅ 1 ⋅
⋅ ⋅ 1
And then there is
You are looking for FillArrays.Eye.
5 Likes
In my early days of Julia I was searching for exactly the same constructor. Somehow the documentation for (I::UniformScaling)(n::Integer) is not visible on docs.julialang.org.
I am not that familiar with Documenter.jl and Git… Does anyone know how we could have it visible in the documentation?
1 Like
The documentation for I can be found in uniformscaling.jl. Please file a pull request adding documentation that I is callable to produce a Diagonal matrix by clicking the “edit” link on that page.
You could also insert a line (LinearAlgebra.I::LinearAlgebra.UniformScaling)(::Integer) after the LinearAlgebra.I line in this index. But I think we would want to cross-reference that from the I documentation anyway.
7 Likes
Just a side-note (EDIT: obsolete, I misread the code, however I left the remaining remark here since it’s valid)
And regarding this:
julia> sizeof(M1), sizeof(M2), sizeof(M3)
(3200, 3200, 8)
M3 stores the diagonal elements in a field called diag. The size you see when you do sizeof(M3) is the address size of your machine (64bit = 8byte), which points to that object:
julia> M3.diag
20-element Array{Int64,1}:
1
1
1
...
...
...
1
1
1
julia> sizeof(M3.diag)
160
As you can see, the elements are in M3.diag with a size of 160, so the overall size is 168 (and not 8, which for 20x “1” would only be possible using some bit juggling and type abstractions).
Whereas M1 and M2 are arrays themselves and sizeof reports the actual size of the whole structure.
1 Like
This has nothing to do with LLVM. fill(1, 20) is a 1d array of 20 ones, so I’m not sure why you think “you basically tell the machine to create a 20x20 matrix, fill it with ones, and then only keep the diagonal elements”. That’s not what this code means in Julia.
Oh, my bad, I read it wrongly, you are perfectly right. @stevengj I removed that part to avoid confusion. Thanks!
Thanks for your response. I did a pull request to improve the docs 
3 Likes
Wait, why did you use Integer for the type here? That’s an abstract type and is what’s causing your problem with this construct. Use Int (or Int64) to fix this.
4 Likes
This (and similar things) happened to me quite a few times already. Sometimes I wish there would some explicit ways to distinguish abstract types (IDE, naming conventions, whatever…)
1 Like
Hm:
julia> M3 = Diagonal(fill(1,20));
julia> M4 = I(20);
julia> sizeof(M3), sizeof(M4)
(8, 8)
julia> sizeof(M3.diag), sizeof(M4.diag)
(160, 20)
Seems like M3 stores integers on the diagonal, while M4 stores booleans on the diagonal.
Yes, as @Tamas_Papp already pointed out in his reply above Linear Algebra tricks?
julia> typeof(M4.diag)
Array{Bool,1}
julia> sizeof(Bool)
1
2 Likes
Ah – @Tamas_Papp’s answer touched on two issues, then. (i) that there is a constructor I(n), and (ii) that this constructor creates a matrix with boolean elements. I overlooked the second issue.
Is there an advantage in using FillArrays.Eye over Diagonal(fill(1,n))? OK, I guess I can install FillArrays.jl and check for myself. But if there is no advantage in my case, I’ll stick with the built-in solution.
OK – I tested it:
julia> M1 = Eye{Int}(5);
julia> M2 = I(5);
julia> M3 = Diagonal(fill(1,5));
julia> sizeof(M1), sizeof(M2), sizeof(M3)
(8, 8, 8)
julia> sizeof(M1.diag), sizeof(M2.diag), sizeof(M3.diag)
(8, 5, 40)
… and then I tested it more:
julia> M1 = rand(-9:9,3,7)
3×7 Array{Int64,2}:
0 -6 1 3 6 -7 9
-4 -8 8 -5 8 -9 9
8 2 -7 -8 0 -5 -5
julia> M2 = [M1 I]
3×10 Array{Int64,2}:
0 -6 1 3 6 -7 9 1 0 0
-4 -8 8 -5 8 -9 9 0 1 0
8 2 -7 -8 0 -5 -5 0 0 1
julia> sizeof(M1)
168
julia> sizeof(M2)
240
julia> M3 = [M1 Eye(3)]
3×10 SparseArrays.SparseMatrixCSC{Float64,Int64} with 22 stored entries:
[2, 1] = -4.0
[3, 1] = 8.0
[1, 2] = -6.0
[2, 2] = -8.0
[3, 2] = 2.0
[1, 3] = 1.0
[2, 3] = 8.0
[3, 3] = -7.0
[1, 4] = 3.0
[2, 4] = -5.0
[3, 4] = -8.0
[1, 5] = 6.0
[2, 5] = 8.0
[1, 6] = -7.0
[2, 6] = -9.0
[3, 6] = -5.0
[1, 7] = 9.0
[2, 7] = 9.0
[3, 7] = -5.0
[1, 8] = 1.0
[2, 9] = 1.0
[3, 10] = 1.0
julia> sizeof(M2), sizeof(M3)
(240, 40)
So – the FillArrays package figured out that M1 has two zeros, so 19 nonzero elements + the three diagonal elements of Eye(3). That is clever.
1 Like
@BLI isn’t this what you were looking for?
julia> Matrix{Int}(I,3,3)
3×3 Array{Int64,2}:
1 0 0
0 1 0
0 0 1
julia> Matrix{Float64}(I,3,3)
3×3 Array{Float64,2}:
1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
julia> Matrix{ComplexF32}(I,3,3)
3×3 Array{Complex{Float32},2}:
1.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 1.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 1.0+0.0im
# etc. etc.
1 Like
You can also use function Base.summarysize which compute the amount of memory, in bytes, used by all unique objects reachable from the argument.
julia> M1 = Eye{Int}(5);
julia> M2 = I(5);
julia> M3 = Diagonal(fill(1,5));
julia> Base.summarysize(M1), Base.summarysize(M2), Base.summarysize(M3)
(8, 53, 88)
and also
julia> M1 = rand(-9:9,3,7);
julia> M2 = [M1 I];
julia> M3 = [M1 Eye(3)];
julia> Base.summarysize(M1), Base.summarysize(M2), Base.summarysize(M3)
(208, 280, 616)
3 Likes