# Custom StaticArray matrix type constructors from sub-arrays of a larger matrix type

**URL:** <https://discourse.julialang.org/t/custom-staticarray-matrix-type-constructors-from-sub-arrays-of-a-larger-matrix-type/83974>\
**Category:** New to Julia\
**Tags:** question, staticarrays\
**Created:** [July 8, 2022, 11:58pm UTC](https://discourse.julialang.org/t/custom-staticarray-matrix-type-constructors-from-sub-arrays-of-a-larger-matrix-type/83974 "2022-07-08T23:58:04Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![mfiano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfiano/32/37565_2.png) [@mfiano](https://discourse.julialang.org/u/mfiano)\
**Post date:** [July 8, 2022, 11:58pm UTC](https://discourse.julialang.org/t/custom-staticarray-matrix-type-constructors-from-sub-arrays-of-a-larger-matrix-type/83974/1 "2022-07-08T23:58:04Z")

</div>

I have a `Matrix2` and `Matrix3` type derived from `StaticArrays.FieldMatrix`:

```julia
(Practice) julia> x = Matrix2(1, 2, 3, 4)
2×2 Matrix2{Int64} with indices SOneTo(2)×SOneTo(2):
 1 3
 2 4

(Practice) julia> y = Matrix3(1, 2, 3, 4, 5, 6, 7, 8, 9)
3×3 Matrix3{Int64} with indices SOneTo(3)×SOneTo(3):
 1 4 7
 2 5 8
 3 6 9

```

I know Julia has some excellent matrix composition functions, and I’m wondering what would be the best way to create a new `Matrix2` outer constructor method such that this is produced:

```julia
(Practice) julia> Matrix2(y)
2×2 Matrix2{Int64} with indices SOneTo(2)×SOneTo(2):
 1 4
 2 5

```

That is, it aligns a window of the 2x2 matrix over the the top left of the 3x3 matrix and extracts those elements, discarding the rest. I’m wondering if this can be done in a general fashion to work with other pairs of square matrices, rather than hard-coding the fields to extract for each constructor method.

I’ve been staring at the docs for a while and not seeing anything too useful here. `reshape` requires the number of elements not to change. `resize!` is for vectors. The only things that looked possibly relevant were `SubArray` types and views into arrays, but I don’t know how to use them with immutable StaticArray types, if it’s even possible.

I need some pointers on what a general solution would be for a fixed set of square matrix types. Ideally I’d like one constructor that does the right thing based on the statically-known square matrix size, rather than writing separate methods with hard-coded fields to extract for each.

Any help would be appreciated. Thanks in advance.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [July 9, 2022, 1:07am UTC](https://discourse.julialang.org/t/custom-staticarray-matrix-type-constructors-from-sub-arrays-of-a-larger-matrix-type/83974/2 "2022-07-09T01:07:40Z")

</div>

Maybe something like this (for it to be efficient, all sizes must be known at compile time):

```julia
julia> function sub_matrix(::Val{N}, m::SMatrix{M,M,T}) where {N,M,T}
           SMatrix{N,N,T}(@view(m[CartesianIndices((1:N,1:N))]))
       end
sub_matrix (generic function with 1 method)

julia> m = SMatrix{3,3,Int,9}(1,2,3,4,5,6,7,8,9)
3×3 SMatrix{3, 3, Int64, 9} with indices SOneTo(3)×SOneTo(3):
 1 4 7
 2 5 8
 3 6 9

julia> sub_matrix(Val(2), m)
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 1 4
 2 5

julia> sub_matrix(Val(1), m)
1×1 SMatrix{1, 1, Int64, 1} with indices SOneTo(1)×SOneTo(1):
 1

julia> sub_matrix(Val(3), m)
3×3 SMatrix{3, 3, Int64, 9} with indices SOneTo(3)×SOneTo(3):
 1 4 7
 2 5 8
 3 6 9

julia> @btime sub_matrix(Val(2), $m)
  6.018 ns (0 allocations: 0 bytes)
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 1 4
 2 5

```

You can also add offsets, to sample some other part of the larger matrix:

```julia
julia> function sub_matrix(::Val{N}, m::SMatrix{M,M,T}; offset=(0,0)) where {N,M,T}
           r1 = 1+offset[1]:N+offset[1]
           r2 = 1+offset[2]:N+offset[2]
           SMatrix{N,N,T}(@view(m[CartesianIndices((r1,r2))]))
       end
sub_matrix (generic function with 1 method)

julia> sub_matrix(Val(2), m)
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 1 4
 2 5

julia> sub_matrix(Val(2), m; offset=(1,1))
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 5 8
 6 9

julia> @btime sub_matrix(Val(2), $m; offset=(1,1))
  14.786 ns (0 allocations: 0 bytes)
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 5 8
 6 9

```

---

<div class="post-metadata">

**Author:** ![mfiano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mfiano/32/37565_2.png) [@mfiano](https://discourse.julialang.org/u/mfiano)\
**Post date:** [July 9, 2022, 3:26am UTC](https://discourse.julialang.org/t/custom-staticarray-matrix-type-constructors-from-sub-arrays-of-a-larger-matrix-type/83974/3 "2022-07-09T03:26:17Z")

</div>

Thank you very much! I’m currently thinking how to extend this to support a `Val` larger than the matrix. Currently `Val(n)` must be \<= `m`, but it would be nice to be able to say copy only the 3x3 portion of a 4x4 matrix into a new 3x3 matrix.

It would also be nice to be able to extend this so that copying a smaller matrix to a larger one leaves the elements along the main diagonal that it doesn’t touch to 1, with the rest 0 (copy a smaller sub-matrix into a larger identity matrix).

These things are going to take me a while to think about. But thank you for the starting point!

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [July 9, 2022, 10:33am UTC](https://discourse.julialang.org/t/custom-staticarray-matrix-type-constructors-from-sub-arrays-of-a-larger-matrix-type/83974/4 "2022-07-09T10:33:48Z")

</div>

You can have a lot of flexibility if you use a mutable static matrix as an intermediate local to the function:

```julia
julia> function newm(::Val{N}, m::SMatrix{M,M,T}) where {N,M,T}
           a = zeros(MMatrix{N,N,T})
           for i in 1:N
               a[i,i] = one(T)
           end
           a[1:M,1:M] .= m
           return SMatrix{N,N,T}(a)
       end
newm (generic function with 1 method)

julia> newm(Val(4), m)
4×4 SMatrix{4, 4, Int64, 16} with indices SOneTo(4)×SOneTo(4):
 1 4 7 0
 2 5 8 0
 3 6 9 0
 0 0 0 1

julia> @btime newm(Val(4), $m)
  11.524 ns (0 allocations: 0 bytes)
4×4 SMatrix{4, 4, Int64, 16} with indices SOneTo(4)×SOneTo(4):
 1 4 7 0
 2 5 8 0
 3 6 9 0
 0 0 0 1

```

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [July 12, 2022, 7:43am UTC](https://discourse.julialang.org/t/custom-staticarray-matrix-type-constructors-from-sub-arrays-of-a-larger-matrix-type/83974/5 "2022-07-12T07:43:31Z")

</div>

Note that it might be somewhat faster to not allocate the intermediate mutable matrix, and construct this by concatenating block matrices (which all produce immutable intermediate types):

```julia
julia> function newm2(::Val{N}, m::SMatrix{M,M,T}) where {N,M,T}
           r = range(0, 0, length=(N-M)*M)
           Z = SMatrix{M,N-M,T,M*(N-M)}(r)
           A = hcat(m, Z)
           B = hcat(Z', Diagonal(SVector{N-M,T}(range(1,1,length=N-M))))
           vcat(A, B)
       end
newm2 (generic function with 1 method)

julia> newm(Val(6), m) == newm2(Val(6), m)
true

julia> @btime newm(Val(6), $m);
  33.430 ns (0 allocations: 0 bytes)

julia> @btime newm2(Val(6), $m);
  19.218 ns (0 allocations: 0 bytes)

```

---

<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:** [July 12, 2022, 8:23am UTC](https://discourse.julialang.org/t/custom-staticarray-matrix-type-constructors-from-sub-arrays-of-a-larger-matrix-type/83974/6 "2022-07-12T08:23:07Z")

</div>

Maybe I’m missing some of the requirements (on a phone), but why not index with `SVectors`?

As in

```julia
ind = SA[1, 2]
y[ind, ind]

```

at least for the simplest case?