Inner product grammar is neater than sum(index) in JuMP modeling, But triggers Warning?

Specific Domains Optimization (Mathematical)
21

This is a “bug” in MutableArithmetics, where we don’t have support for SubArray. The biggest problem in Julia is that there are no formal interfaces, so things work until they don’t, and then you hit a MethodError.

I’ve opened an issue: similar_array_type with SubArray · Issue #316 · jump-dev/MutableArithmetics.jl · GitHub

2 Likes
22

I was insane💀, since I could have written:

sum(p .* M_vector)

or in my “another example”

sum(p .* eachslice(M_Array; dims = 3))

which were concise, although still a slight bit of inconvenient compared to

ip(p, M_vector)

But after all, the MethodError in MutableArithmetics is revealed💀.

23

I’m fixing MA: Fix similar_array_type for SubArray by odow · Pull Request #317 · jump-dev/MutableArithmetics.jl · GitHub

2 Likes
24

When I run ip(p, M_vector), I get

ERROR: MethodError: no method matching zero(::Type{Matrix{AffExpr}})
The function `zero` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  zero(::Type{Union{}}, Any...)
   @ Base number.jl:310
  zero(::Type{Dates.Time})
   @ Dates ~/.julia/juliaup/julia-1.11.2+0.x64.linux.gnu/share/julia/stdlib/v1.11/Dates/src/types.jl:460
  zero(::Type{Pkg.Resolve.VersionWeight})
   @ Pkg ~/.julia/juliaup/julia-1.11.2+0.x64.linux.gnu/share/julia/stdlib/v1.11/Pkg/src/Resolve/versionweights.jl:15
  ...

Stacktrace:
 [1] neutral_element(::typeof(+), T::Type)
   @ MutableArithmetics ~/.julia/dev/MutableArithmetics/src/reduce.jl:25
 [2] fused_map_reduce(::typeof(add_dot), ::Vector{VariableRef}, ::Vector{Matrix{Float64}})
   @ MutableArithmetics ~/.julia/dev/MutableArithmetics/src/reduce.jl:44
 [3] operate(::typeof(dot), x::Vector{VariableRef}, y::Vector{Matrix{Float64}})
   @ MutableArithmetics ~/.julia/dev/MutableArithmetics/src/implementations/LinearAlgebra.jl:508
 [4] dot(x::Vector{VariableRef}, y::Vector{Matrix{Float64}})
   @ MutableArithmetics ~/.julia/dev/MutableArithmetics/src/dispatch.jl:58

I have opened an issue dot product fails when the result is a matrix · Issue #318 · jump-dev/MutableArithmetics.jl · GitHub

25

Does it mean that after this revision, I could write

transpose(p) * [rand(2, 2) for _ in 1:length(p)] # p::Vector{JuMP.VariableRef}

without seeing an red ERROR messag?
This transpose(p) * v Grammar is also likeable, a bit better than sum(p .* v)

26

Yes, transpose(p) * v redirects to LinearAlgebra._dot_nonrecursive which is then caught by MutableArithmetics because eltype(p) is a JuMP variable and so it’s then redirected to MutableArithmetics.fused_map_reduce(MutableArithmetics.add_mul, p, v) which is working now with Fix fused_map_reduce when result is a matrix by blegat · Pull Request #319 · jump-dev/MutableArithmetics.jl · GitHub

1 Like
27

We just released MutableArithmetics v1.6.1 with the fix

28

I reviewed some math books and concluded that ip (inner product) is a false name which conveyed a misconception. A correction could be

function sm(p, x) return sum(p .* x) end

where s stands for vector sum and m stands for scalar multiplication.

The sm should be a more fundamental operation, which makes sense in, e.g.,
integration theory (p represents a probability measure while x represents a measurable function).

And some related useful formulas are:

  1. @assert sm(A, B) == tr(transpose(A) * B)
  2. @assert transpose(z) * A * w == tr(A * w * transpose(z))
  3. @assert tr(A * B) == tr(B * A) as long as size(A) == size(transpose(B))
  4. @assert sm(c .* x, y) == sm(c, diag(x * y'))

where, tr and diag are from module LinearAlgebra.
From the 4th point we see that sm is indeed more versatile that tr. Therefore I think it’s obfuscating to define an inner product using tr (as it was done in some textbooks).

29

Given the discussion here. I think we are motivated to fathom whether there is a standard (or at least recommended) way to write JuMP code that attains the optimal performance. @odow

Some very common cases include

import JuMP
using LinearAlgebra
model = JuMP.Model();
JuMP.@variable(model, x[1:3]);
JuMP.@variable(model, X[1:3, 1:3]);

JuMP.@expression(model, lin_scalar_1, rand(3)' * x)
JuMP.@expression(model, lin_scalar_2, dot(rand(3), x))

JuMP.@expression(model, quad_scalar_1, x' * rand(3, 3) * x)
JuMP.@expression(model, quad_scalar_2, dot(x, rand(3, 3), x))

JuMP.@expression(model, matrix_lin_scalar_1, rand(3)' * X * rand(3))
JuMP.@expression(model, matrix_lin_scalar_2, dot(rand(3), X, rand(3)))

JuMP.@expression(model, mat_lin_scalar_1, sum(rand(3, 3) .* X))
JuMP.@expression(model, mat_lin_scalar_2, dot(rand(3, 3), X))

It appears that they will have different dispatches in MutableArithmetics, thus possibly different performances.

For a reference, notice the word Unlike in the following docstring

help?> *

  *(A, B::AbstractMatrix, C)
  A * B * C * D

  If the last factor is a vector, or the first a transposed vector, then it is efficient to deal with      
  these first. In particular x' * B * y means (x' * B) * y for an ordinary column-major B::Matrix. Unlike  
  dot(x, B, y), this allocates an intermediate array.