# Linear Algebra tricks?

**URL:** <https://discourse.julialang.org/t/linear-algebra-tricks/38618>\
**Category:** General Usage\
**Created:** [May 2, 2020, 5:15pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618 "2020-05-02T17:15:05Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 2, 2020, 5:15pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/1 "2020-05-02T17:15:05Z")

</div>

Two questions related to Linear Algebra:

1. I need an identity matrix in some computations. Two ways to achieve this:

```julia
julia> Matrix{Integer}(I,3,3)
3×3 Array{Integer,2}:
 true 0 0
    0 true 0
    0 0 true

```

or

```julia
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
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?]

---

<div class="post-metadata">

**Author:** ![jbrea](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jbrea/32/3879_2.png) [@jbrea](https://discourse.julialang.org/u/jbrea)\
**Post date:** [May 2, 2020, 5:19pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/2 "2020-05-02T17:19:19Z")

</div>

> [@BLI](#):
>
> _Question 1_ : what is the preferred method?

Are you aware of

```julia
help?> I

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

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 2, 2020, 5:44pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/3 "2020-05-02T17:44:36Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [May 2, 2020, 5:45pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/4 "2020-05-02T17:45:03Z")

</div>

You can use Diagonal?

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 2, 2020, 5:48pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/5 "2020-05-02T17:48:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 2, 2020, 5:57pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/6 "2020-05-02T17:57:14Z")

</div>

Wow!

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

```

---

<div class="post-metadata">

**Author:** ![CodeLenz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/codelenz/32/3419_2.png) [@CodeLenz](https://discourse.julialang.org/u/CodeLenz)\
**Post date:** [May 2, 2020, 6:09pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/7 "2020-05-02T18:09:04Z")

</div>

It is possible to use

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

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 2, 2020, 6:21pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/8 "2020-05-02T18:21:50Z")

</div>

Superb! So both of my questions are answered, I guess:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [May 2, 2020, 6:25pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/9 "2020-05-02T18:25:41Z")

</div>

> [@BLI](#):
>
> Yes, but there are cases when one needs to specify the size.

There is an `I(n)` constructor, eg

```julia
julia> I(3)
3×3 Diagonal{Bool,Array{Bool,1}}:
 1 ⋅ ⋅
 ⋅ 1 ⋅
 ⋅ ⋅ 1

```

And then there is

> **[GitHub - JuliaArrays/FillArrays.jl: Julia package for lazily representing...](https://github.com/JuliaArrays/FillArrays.jl)**
>
> Julia package for lazily representing matrices filled with a single entry - GitHub - JuliaArrays/FillArrays.jl: Julia package for lazily representing matrices filled with a single entry

You are looking for `FillArrays.Eye`.

---

<div class="post-metadata">

**Author:** ![roble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roble/32/12303_2.png) [@roble](https://discourse.julialang.org/u/roble)\
**Post date:** [May 2, 2020, 9:29pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/10 "2020-05-02T21:29:17Z")

</div>

In my early days of Julia I was searching for exactly the same constructor. Somehow the documentation for [`(I::UniformScaling)(n::Integer)`](https://github.com/JuliaLang/julia/blob/381693d3dfc9b7072707f6d544f82f6637fc5e7c/stdlib/LinearAlgebra/src/uniformscaling.jl#L74) is not visible on [docs.julialang.org](http://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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 2, 2020, 9:46pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/11 "2020-05-02T21:46:49Z")

</div>

> [@roble](#):
>
> Does anyone know how we could have it visible in the documentation?

The documentation for `I` can be [found in `uniformscaling.jl`](https://github.com/JuliaLang/julia/blob/7be230bb6178ed74297a92d2e50ef8d8d801351d/stdlib/LinearAlgebra/src/uniformscaling.jl#L33-L48). Please file a pull request adding documentation that `I` is callable to produce a `Diagonal` matrix by [clicking the “edit” link](https://help.github.com/en/github/managing-files-in-a-repository/editing-files-in-your-repository) on that page.

You could also insert a line `(LinearAlgebra.I::LinearAlgebra.UniformScaling)(::Integer)` after the `LinearAlgebra.I` line in [this index](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/docs/src/index.md). But I think we would want to cross-reference that from the `I` documentation anyway.

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [May 2, 2020, 10:39pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/12 "2020-05-02T22:39:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 2, 2020, 10:43pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/13 "2020-05-02T22:43:43Z")

</div>

> [@tamasgal](#):
>
> To be honest, I first thought `Diagonal(fill(1, 20))` will in fact create the matrix first, but luckily LLVM is smarter than me 😉

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.

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [May 2, 2020, 10:44pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/14 "2020-05-02T22:44:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![roble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roble/32/12303_2.png) [@roble](https://discourse.julialang.org/u/roble)\
**Post date:** [May 2, 2020, 11:08pm UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/15 "2020-05-02T23:08:30Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 3, 2020, 12:20am UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/16 "2020-05-03T00:20:03Z")

</div>

> [@BLI](#):
>
> ```julia
> Matrix{Integer}(I,3,3)
> 
> ```

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.

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [May 3, 2020, 9:55am UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/17 "2020-05-03T09:55:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 3, 2020, 10:43am UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/18 "2020-05-03T10:43:26Z")

</div>

Hm:

```julia
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.

---

<div class="post-metadata">

**Author:** ![tamasgal](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamasgal/32/27946_2.png) [@tamasgal](https://discourse.julialang.org/u/tamasgal)\
**Post date:** [May 3, 2020, 10:45am UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/19 "2020-05-03T10:45:19Z")

</div>

Yes, as @Tamas_Papp already pointed out in his reply above [Linear Algebra tricks? - #9 by Tamas\_Papp](https://discourse.julialang.org/t/linear-algebra-tricks/38618/9)

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

julia> sizeof(Bool)
1

```

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [May 3, 2020, 10:52am UTC](https://discourse.julialang.org/t/linear-algebra-tricks/38618/20 "2020-05-03T10:52:46Z")

</div>

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.

[Next page](https://discourse.julialang.org/t/linear-algebra-tricks/38618.md?page=2)