# Linear Algebra tricks?

Two questions related to Linear Algebra:

1. I need an identity matrix in some computations. Two ways to achieve this:
``````julia> Matrix{Integer}(I,3,3)
3×3 Array{Integer,2}:
true     0     0
0  true     0
0     0  true
``````

or

``````julia> diagm(fill(1,3))
3×3 Array{Int64,2}:
1  0  0
0  1  0
0  0  1
``````

I suspect the first approach is more memory efficient (?), but it looks ugly with the boolean interpretation of `1`. It gets even worse if I create an identity matrix over the `Rational` type… (`true//true`)

Question 1: what is the preferred method?

1. I’m doing some column permutations on a matrix `A`, collected in vector `p`. Thus, the resulting matrix is `A[:,cp]`. There is an equivalent permutation matrix `P` so that `A*P = A[:,cp]`. (Essentially, `P = diagm(fill(1,n))[:,p]` where `n=size(A,2)`).
I now need to do row permutations on another matrix, say `P*B*`. So I need to convert the column permutation `cp` vector to a row permutation vector `rp` so that `P*B = B[rp,:]`.

Question 2: Is there a built-in Julia function for doing this conversion from `cp` to `rp`?

[So far, I’ve done it as follows:

``````julia> A = rand(-9:9,3,3)
3×3 Array{Int64,2}:
0  4  -3
5  2  -5
-6  7   0

julia> cp = [3,1,2]
3-element Array{Int64,1}:
3
1
2

julia> A1 = A[:,cp]
3×3 Array{Int64,2}:
-3   0  4
-5   5  2
0  -6  7

julia> P = diagm(fill(1,3))[:,cp]
3×3 Array{Int64,2}:
0  1  0
0  0  1
1  0  0

julia> A*P
3×3 Array{Int64,2}:
-3   0  4
-5   5  2
0  -6  7

julia> P*A
3×3 Array{Int64,2}:
5  2  -5
-6  7   0
0  4  -3

julia> rp = P*(1:3)
3-element Array{Int64,1}:
2
3
1

julia> A[rp,:]
3×3 Array{Int64,2}:
5  2  -5
-6  7   0
0  4  -3
``````

So… `rp = diagm(fill(1,n))[:,cp]*(1:n)`. But this is perhaps not very efficient? Is there a better way?]

1 Like

Are you aware of

``````help?> I

An object of type UniformScaling, representing an identity matrix of any size.
``````
5 Likes

Yes, but there are cases when one needs to specify the size.

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.