# Extend for vector-valued matrix multiplication?

**URL:** <https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286>\
**Category:** General Usage\
**Tags:** question, array\
**Created:** [March 21, 2020, 3:31am UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286 "2020-03-21T03:31:58Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [March 21, 2020, 3:31am UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/1 "2020-03-21T03:31:58Z")

</div>

Consider a `UV::Vector{SVector{2}}` or alternatively `(U,V)::Tuple{Vector,Vector}`.

I would like to compute the same matrix multiplication formula for `U` and `V` but would prefer storing the information in the `Vector{SVector{2}}` format (technically, it’s not an `SVector` but my own type).

```nohighlight
U = U - M*U
V = V - M*V

```

My preference would be for a single formula `UV = UV - M*UV` where `UV::Vector{SVector{2}}` is a vector of static vectors (or my own similar static type actually).

The matrix `M` is a sparse matrix with only `Real` coefficients. In order for `M*UV` to make sense, the matrix multiplication needs to broadcast each operation onto the `SVector{2}` in order to do the two matrix multiplications simultaneously in a single equation instead of two equations.

The point is that it is much better for me to have `UV` stored in memory, than to have `U,V` separate. It would save me some allocations elsewhere if I could do the matrix multiplication like this.

How would someone go about extending the matrix multiplication for this purpose?

---

<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:** [March 21, 2020, 5:45am UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/2 "2020-03-21T05:45:38Z")

</div>

As a sketch you could e.g. make a new type `UV` and then make a new method `*(M, uv::UV)` for `*`.

Basically provided there is at least one of your own types in there, you can redefine any operations you want. What you shouldn’t do is redefine operations on only base types (this is “type piracy”).

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [March 21, 2020, 2:41pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/3 "2020-03-21T14:41:42Z")

</div>

As stated in my post, I’m using my own types, but for familiarity I used common types in my post.

---

<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:** [March 21, 2020, 3:13pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/4 "2020-03-21T15:13:24Z")

</div>

```julia
julia> using StaticArrays

julia> Base.:*(M::AbstractMatrix, v::Vector{<:SVector}) = [M*x for x in v]

julia> UV = [SA[1, 2], SA[3, 4]]
2-element Array{SArray{Tuple{2},Int64,1,2},1}:
[1, 2]
[3, 4]

julia> M * UV
2-element Array{Array{Int64,1},1}:
[5, 11]
[11, 25]

```

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [March 21, 2020, 3:16pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/5 "2020-03-21T15:16:02Z")

</div>

You are misunderstanding the layout of my problem, it’s more like this:

```nohighlight
UV = [rand(SVector{2,Float64}) for k ∈ 1:1000]
M = rand(1000,1000)

```

---

<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:** [March 21, 2020, 3:16pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/6 "2020-03-21T15:16:26Z")

</div>

But it’s still probably better to define your own type

```julia
struct UV{T}
    v::Vector{SVector{2,T}}
end

```

or better as an `SVector` of `SVector`s perhaps.

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [March 21, 2020, 3:21pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/7 "2020-03-21T15:21:19Z")

</div>

What I need is something like this

```julia
[M*getindex.(UV,i) for i in 1:2]

```

But the result needs to be a `Vector{SVector{2}}` layout, I’d like to do this without re-allocating.

---

<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:** [March 21, 2020, 3:59pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/8 "2020-03-21T15:59:12Z")

</div>

> [@chakravala](#):
>
> You are misunderstanding the layout of my problem, it’s more like this:
> 
> ```julia
> UV = [rand(SVector{2,Float64}) for k ∈ 1:1000]
> M = rand(1000,1000)
> 
> ```

I’m misunderstanding because you didn’t specify the problem well. In this case you can’t multiply that matrix by each element of `UV`?

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [March 21, 2020, 4:18pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/9 "2020-03-21T16:18:47Z")

</div>

No, each element of `UV` is of `length` 2 but `M` is e.g. 1000\times1000, so multiplying `M` by each element of `UV` is not the goal, since that would be a dimension mismatch.

---

<div class="post-metadata">

**Author:** ![bcmichael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bcmichael/32/10003_2.png) [@bcmichael](https://discourse.julialang.org/u/bcmichael)\
**Post date:** [March 21, 2020, 4:38pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/10 "2020-03-21T16:38:42Z")

</div>

I think your desire to not reallocate is going to be a problem for a couple of reasons. First of all it’s not really possible to do a matrix-vector multiplication and store the results in the same vector you’re multiplying by; you really need to store the result somewhere different. Second of all SVectors are immutable, so trying to modify them will be a problem.

Julia isn’t doing anything super special like calling into external libraries to do linear algebra with sparse matrices. The matrix multiplication routines are just pretty simple functions defined in [https://github.com/JuliaLang/julia/blob/master/stdlib/SparseArrays/src/linalg.jl](https://github.com/JuliaLang/julia/blob/master/stdlib/SparseArrays/src/linalg.jl)

I think the basic approach for your function is going to have to entail a separate array to store the results in then make new SVectors with the results to stick back into UV. As long as you can preallocate that temporary results array I don’t think that will really entail actual allocation as the compiler should be pretty good about handling static vectors. Does adapting the first function definition in that file to something like the following do what you want?

```julia
using StaticArrays
using SparseArrays
import SparseArrays.getcolptr

struct UV{T}
    v::Vector{SVector{2,T}}
end

function my_multiplication_operation(C::Array{T,2}, A::SparseMatrixCSC, B::UV{T}) where T
    size(A, 2) == size(B.v, 1) || throw(DimensionMismatch())
    size(A, 1) == size(A, 2) || throw(DimensionMismatch())
    size(A, 1) == size(C, 1) || throw(DimensionMismatch())
    size(C, 2) == 2 || throw(DimensionMismatch())
    nzv = nonzeros(A)
    rv = rowvals(A)
    fill!(C, zero(eltype(C)))
    for k = 1:size(C, 2)
        @inbounds for col = 1:size(A, 2)
            αxj = B.v[col][k]
            for j = getcolptr(A)[col]:(getcolptr(A)[col + 1] - 1)
                C[rv[j], k] += nzv[j]*αxj
            end
        end
    end
    for i = 1:size(C,1)
        B.v[i] = SVector{2,T}(B.v[i][1]-C[i,1], B.v[i][2]-C[i,2])
    end
end

a = sprand(6,6,0.3)
b = rand(6,2)
uv = UV([SVector{2}(b[n,1], b[n,2]) for n = 1:6])
c = Array{Float64,2}(undef, (6,2))

d = b-a*b
my_multiplication_operation(c,a,uv)
all((uv.v[n][1] == d[n,1] && uv.v[n][2] == d[n,2] for n in 1:6))

```

---

<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:** [March 21, 2020, 4:52pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/11 "2020-03-21T16:52:17Z")

</div>

Right, so please specify precisely the actual goal?

---

<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:** [March 21, 2020, 4:54pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/12 "2020-03-21T16:54:24Z")

</div>

@chakravala Do you just want something like this?

```julia
julia> using StaticArrays

julia> input = rand(SVector{2,Float64}, 1000);

julia> output = similar(input); # allocation one time

julia> M = SA[1 2; 3 4];

julia> output .= (M * x for x in input); # no new allocation

```

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [March 21, 2020, 5:25pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/13 "2020-03-21T17:25:55Z")

</div>

@dpsanders no that’s not what I’m trying to do, `M` is not the same dimension as `x`.

@bcmichael thanks, that’s similar to what I want, here is a modification of it with the result I want

```nohighlight
function my_multiplication_operation(A::SparseMatrixCSC, B::UV{T}) where T
    C = zeros(SVector{2,T},length(B.v))
    nzv = nonzeros(A)
    rv = rowvals(A)
    @inbounds for col = 1:length(B.v)
        αxj = B.v[col]
        for j = getcolptr(A)[col]:(getcolptr(A)[col + 1] - 1)
            C[rv[j]] += nzv[j]*αxj
        end
    end
    return C
end

```

Let me know if you know have any further suggestion to improve this.

```nohighlight
julia> my_multiplication_operation(a,uv)
6-element Array{SArray{Tuple{2},Float64,1,2},1}:
 [0.469116188764566, -0.06476447346330405]  
 [0.0, 0.0]                                 
 [0.6690022529288051, 0.5049880060027236]   
 [0.0, 0.0]                                 
 [-0.41935221191364325, -0.5347303406796164]
 [0.1903883495009893, 0.7041576907472157]   

```

Is there any detailed explanation of `getcolptr` available?

---

<div class="post-metadata">

**Author:** ![bcmichael](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bcmichael/32/10003_2.png) [@bcmichael](https://discourse.julialang.org/u/bcmichael)\
**Post date:** [March 21, 2020, 6:30pm UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/14 "2020-03-21T18:30:59Z")

</div>

The definition of the fields in SparseMatrixCSC looks like

```julia
struct SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrixCSC{Tv,Ti}
    m::Int # Number of rows
    n::Int # Number of columns
    colptr::Vector{Ti} # Column i is in colptr[i]:(colptr[i+1]-1)
    rowval::Vector{Ti} # Row indices of stored values
    nzval::Vector{Tv} # Stored values, typically nonzeros

```

Looks like in this case getcolptr just returns the colptr field, but it has another method to deal with views of sparse matrices.

---

<div class="post-metadata">

**Author:** ![chakravala](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chakravala/32/6832_2.png) [@chakravala](https://discourse.julialang.org/u/chakravala)\
**Post date:** [March 25, 2020, 7:37am UTC](https://discourse.julialang.org/t/extend-for-vector-valued-matrix-multiplication/36286/15 "2020-03-25T07:37:14Z")

</div>

In the end, all I needed was this short definition

> <https://github.com/chakravala/Grassmann.jl/blob/e242eba2b288371ce7eb756d543c653d1c1ef7b6/src/composite.jl#L385-L386>