Linear Algebra tricks?

General Usage
1

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
2

Are you aware of

help?> I

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

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

4

You can use Diagonal?

3 Likes
5

OK – for a vector v, diagm(v) and Diagonal(v) give the same result? Or does Diagonal(v) avoid storing off-diagonal elements?

6

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)
7

It is possible to use

A = I(2)
2×2 Diagonal{Bool,Array{Bool,1}}:
1 ⋅
⋅ 1

4 Likes
8

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
9

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
10

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
11

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
12

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
13

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.

14

Oh, my bad, I read it wrongly, you are perfectly right. @stevengj I removed that part to avoid confusion. Thanks!

15

Thanks for your response. I did a pull request to improve the docs :white_check_mark:

3 Likes
16

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
17

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
18

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.

19

Yes, as @Tamas_Papp already pointed out in his reply above Linear Algebra tricks? - #9 by Tamas_Papp

julia> typeof(M4.diag)
Array{Bool,1}

julia> sizeof(Bool)
1
2 Likes
20

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.