Allocation of an identity block into a larger matrix

user_9576551 0
Viewed 133

In matlab, assuming I want to build a matrix A containing an identity block, I may use

A = zeros(4,4);
A(1:2,1:2) = eye(2);

in julia, I know of the UniformScaling / I operator. However I does not seem to be usable in the same fashion, and everything I read says that it supersedeed the eye operator.

What is the julianic way to do this ?

EDIT: I am currently using some Diagonal(ones(...)) to actually create the diagonal block.

With respect to the first answer, it seems relatively comparable to the Matrix(I...) syntax (twice as slower max), but much better in terms of memory usage for large matrices (about a hundred times less).

See this code for test purposes

using LinearAlgebra
using BenchmarkTools


A = zeros(10,10)
@btime A[1:8,1:8] = Matrix(I,8,8)

B = zeros(10,10)
@btime B[1:8,1:8] = Diagonal(ones(8))

A == B

A = zeros(100,100)
@btime A[1:80,1:80] = Matrix(I,80,80)

B = zeros(100,100)
@btime B[1:80,1:80] = Diagonal(ones(80))

A == B

A = zeros(1000,1000)
@btime A[1:800,1:800] = Matrix(I,800,800)

B = zeros(1000,1000)
@btime B[1:800,1:800] = Diagonal(ones(800))

A == B

which returns

  155.858 ns (1 allocation: 160 bytes)
  176.879 ns (2 allocations: 160 bytes)
  6.680 μs (1 allocation: 6.44 KiB)
  10.799 μs (2 allocations: 752 bytes)
  617.099 μs (2 allocations: 625.14 KiB)
  1.007 ms (2 allocations: 6.39 KiB)
true
julia
2 Answers

Here's a generic way that hardly allocates at all:

julia> using LinearAlgebra, BenchmarkTools

julia> function func1(z, a, b)
           z[diagind(z)[a:b]] .= one(eltype(z))
           return z   end
func1 (generic function with 1 method)

julia> test = zeros(1_000, 1_000);

julia> @btime func1($test, 1, 2);
  68.680 ns (2 allocations: 80 bytes)

julia> @btime func1($test, 1, 500);
  441.919 ns (2 allocations: 80 bytes)

julia> @btime func1($test, 1, 1000);
  2.444 μs (2 allocations: 80 bytes)

Edited to add: in contrast to many other languages, the simplest and most efficient solution in Julia is often a loop:

julia> function func2(z, a, b)
           for i ∈ a:b
               z[i, i] = one(eltype(z))
           end
           return z
       end
func2 (generic function with 1 method)
julia> @btime func2($test, 1, 2);
  2.000 ns (0 allocations: 0 bytes)

julia> @btime func2($test, 1, 500);
  265.337 ns (0 allocations: 0 bytes)

julia> @btime func2($test, 1, 1000);
  2.200 μs (0 allocations: 0 bytes)
user_24998920

I'm new to Julia myself so feel free to correct me:

julia> z = zeros(Int, 4, 4)
4×4 Matrix{Int64}:
 0  0  0  0
 0  0  0  0
 0  0  0  0
 0  0  0  0

julia> using LinearAlgebra

julia> z[2:3, 2:3] = Matrix{Int}(I, 2, 2)
2×2 Matrix{Int64}:
 1  0
 0  1

julia> z
4×4 Matrix{Int64}:
 0  0  0  0
 0  1  0  0
 0  0  1  0
 0  0  0  0
user_129982050
Related