Matrix-vector multiplication is a bit off from "equivalent" [Strided]Matrix-Strided]Matrix multiplication (and I guess MATLAB)

General Usage
1

I noticed an odd thing when porting code from MATLAB, that this line makes a difference:

 julia> gamma_gls_hat = gamma_gls_hat[:, :]
12×1 Matrix{Float64}

gamma_gls_hat as I calculated before this line was a 12-element Vector{Float64} and MATLAB doesn’t have that concept (at least originally, in MATLAB everything was a matrix, not sure if that changed, allowing vectors now).

So, as this has sometimes been important while translating (Julia is more strict), I added that line after my calculation, and was now checking if I could take it out.

I get the same 12 numbers in the matrix as in the vector (i.e. .== shows that, ==, or === doesn’t or shouldn’t), while none of the values in subsequent calculation are the same:

julia> (tmp*gamma_gls_hat) .== (tmp*gamma_gls_hat[:, :])

julia> (tmp*gamma_gls_hat) - (tmp*gamma_gls_hat[:, :])
5×1 Matrix{Float64}:
 -4.336808689942018e-19
  8.673617379884035e-19
 -2.168404344971009e-19
  2.168404344971009e-19
  3.959879028413854e-20

julia> (tmp*gamma_gls_hat) ≈ (tmp*gamma_gls_hat[:, :])  # still:
true

julia> tmp
5×12 Matrix{Float64}:

julia> VERSION
v"1.6.0-DEV.1177"

I realize it’s only a bit off (and probably would be in all cases), caused by some optimization:

@edit (tmp*gamma_gls_hat[:, :])  # gets me to line 158:

# optimization for dispatching to BLAS, e.g. *(::Matrix{Float32}, ::Matrix{Float64})
# but avoiding the case *(::Matrix{<:BlasComplex}, ::Matrix{<:BlasReal})
# which is better handled by reinterpreting rather than promotion
function (*)(A::StridedMatrix{<:BlasReal}, B::StridedMatrix{<:BlasFloat})

vs.

@edit (tmp*gamma_gls_hat)  # gets me to line 44:

function (*)(A::StridedMatrix{T}, x::StridedVector{S}) where {T<:BlasFloat,S<:Real}

I guess I shouldn’t strive for bit-identical, at least if I wanted, compared to MATLAB, I would need this extra line.

2

The short answer is don’t worry about it and use the vector. The matrix matrix and matrix vector algorithms are different leading to slightly different results, but if anything the matrix vector should be a bit more accurate on average, and are also probably faster.

3 Likes