Odd behavior for matrix multiplication with custom scalar

I was just playing around with the standard trick of computing shortest paths through analogy to matrix exponentiation, and I came across some unexpected behavior. (Unexpected to me, that is.)

The idea is to replace product and sum with sum and minimum, and so I just wrapped the distances in my own distance scalar:

import Base: zero, +, *

struct Dist{T}

zero(x::Dist) = Dist(zero(x.d))
+(x::Dist, y::Dist) = Dist(min(x.d, y.d))
*(x::Dist, y::Dist) = Dist(x.d + y.d)

Now, let’s try to use it:

n = 3
G = Dist.(rand(n, n))
println([x.d for x in G^n])

This works well enough, but once you set n to 4 or higher, all entries become zero. I thought it might be some kind of strange overflow issue or something (sounded unlikely, but still…) so I replaced G^n with G * G, and the same happens. For n = 4, the product zeros out the matrix, but for n = 3, it does not.

I guess this must be a cutoff for some other product algorithm or something (haven’t studied the relevant code yet), but … shouldn’t this still work? Or am I going about it all wrong? Or missing something else obvious?-D

(Note that this is not an important problem for me to solve, subject-matter-wise; I’m just curious about the language aspect, here, and why the custom operations mess with the matrix product.)

For 2x2 and 3x3 specialized methods exist https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/matmul.jl#L752.

For larger matrices a generic one is used https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/matmul.jl#L583.

Right. And this general one fails for custom scalars, or?

1 Like

It’s not supposed to…

help?> zero

  Get the >>>>additive identity<<<<< element for the type of x (x can also specify the type itself).
julia> D = Dist(rand())

julia> D + zero(D)

Heh, I just got the error that zero wasn’t implemented, and added it without considering the role it played. Should have been inf, of course. Thanks!

1 Like