Extend for vector-valued matrix multiplication?

General Usage
,
1

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

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?

1 Like
2

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

2 Likes
3

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

1 Like
4 
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]
1 Like
5

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

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

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

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

or better as an SVector of SVectors perhaps.

1 Like
7

What I need is something like this

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

8

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?

1 Like
9

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.

10

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

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?

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))
1 Like
11

Right, so please specify precisely the actual goal?

12

@chakravala Do you just want something like this?

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
13

@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

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.

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?

14

The definition of the fields in SparseMatrixCSC looks like

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.

15

In the end, all I needed was this short definition