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

1 Like