Counter-intuitive difference in matrix-vector multiplication when the matrix is transposed

Specific Domains Numerics
1

I found that matrix-vector multiplication in Julia uses a different algorithm when the matrix is wrapped by Transpose or Adjoint, compared to when it is a regular Matrix. As a result, the multiplication can sometimes produces slightly different results. You can see this by running the following small script:

A = rand(Float64, (10, 20))
A_transpose = copy(A')
x = rand(size(A, 2))

y1 = A*x
y2 = A_transpose'x

display(hcat(y1, y2, y1 .- y2))

When running this script multiple times, I occasionally get outputs like these:

$ julia test.jl 
10Ă—3 Matrix{Float64}:
 4.04039  4.04039  0.0
 5.88688  5.88688  0.0
 6.64932  6.64932  0.0
 4.84645  4.84645  0.0
 4.87441  4.87441  0.0
 6.07819  6.07819  0.0
 6.5338   6.5338   0.0
 4.84789  4.84789  0.0
 4.56281  4.56281  0.0
 5.19563  5.19563  0.0

$ julia test.jl 
10Ă—3 Matrix{Float64}:
 3.62809  3.62809   0.0
 4.65584  4.65584   8.88178e-16
 3.10216  3.10216   0.0
 3.51676  3.51676  -4.44089e-16
 3.72409  3.72409   0.0
 2.98007  2.98007  -4.44089e-16
 2.43762  2.43762   0.0
 3.16095  3.16095   0.0
 3.87464  3.87464   8.88178e-16
 4.62357  4.62357  -8.88178e-16

Is this difference intentionally introduced, or is it a bug?

2 
import LinearAlgebra
A = rand(10, 20)
x = rand(20)

y1 = A * x
y2 = transpose(A') * x
y3 = (transpose(A))' * x

@assert LinearAlgebra.norm(y1 .- y2) == 0
@assert LinearAlgebra.norm(y3 .- y2) == 0

This hardly error.

3

Neither. That is just floating point rounding errors.
As a result of their finite precision, floating point operations are not associative, e.g.

a=1
b=1e-20
@assert (b+a)-a != b+(a-a)

This means that results naturally depend on the exact order of operation you perform. Results are generally expected to change across algorithms. However for results shouldn’t change much if the algorithms are well designed. What you see here is on the level of the accuracy of Float64. So basically a difference of 1 unit in last place.

If your algorithm/application is sensitive to these tiny deviations, then you should consider

  1. being very very careful with every operation
  2. Use bigger/more precise floats (if suitable)
  3. Use some exact arithmetic (maybe fixed point numbers or rationals)
4 Likes
4

An the reason things are done in a different order when using a transposed matrix is probably due to questions of memory-locality. I.e. it’s usually good for performance to access memory in sequential order.

2 Likes
5

Thanks everyone for the replies.

I already knew the difference came from floating-point arithmetic, but I wanted to check whether it was expected behavior or something unintended. From your answers, I now understand it’s neither an intentional spec nor a bug, and that I should handle it by reordering calculations or making the system more stable.

@abraemer said:

If your algorithm/application is sensitive to these tiny deviations…

Your guess was spot on. When I posted the first message, I was struggling to port a Python implementation to Julia and was trying to figure out why the results weren’t matching. After digging into it, I found the issue was caused by the large condition number of the design matrix and the Jacobian.