# Linear Algebra tricks?

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