# Matrix multiplication and vector multiplication give different results!

**URL:** https://discourse.julialang.org/t/matrix-multiplication-and-vector-multiplication-give-different-results/130177
**Category:** Numerics
**Created:** [June 24, 2025, 4:27pm UTC](https://discourse.julialang.org/t/matrix-multiplication-and-vector-multiplication-give-different-results/130177 "2025-06-24T16:27:32Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)
#### Post date: [June 24, 2025, 4:27pm UTC](https://discourse.julialang.org/t/matrix-multiplication-and-vector-multiplication-give-different-results/130177/1 "2025-06-24T16:27:32Z")

</div>

I found that (slightly) different results coming from matrix-matrix vs. matrix-vector multiplication!

```julia
julia> m1 = randn(200, 300);

julia> m2 = randn(300, 400);

julia> m12 = m1 * m2; # matrix-matrix multiplication

julia> m12col400 = m1 * m2[:, 400]; # matrix-vector multiplication

julia> m12[:, 400] == m12col400
false

julia> m12[:, 400] .- m12col400
200-element Vector{Float64}:
  0.0
  0.0
 -1.4210854715202004e-14
  1.7416623698807143e-15
 -7.521760991835436e-15
  1.2434497875801753e-14
  4.440892098500626e-15
  0.0
  2.6645352591003757e-15
  5.329070518200751e-15
  3.552713678800501e-15
  1.7763568394002505e-15
  3.552713678800501e-15
  ⋮
  3.552713678800501e-15
 -7.105427357601002e-15
  1.3322676295501878e-15
 -3.552713678800501e-15
  3.1086244689504383e-15
 -1.7763568394002505e-15
 -3.552713678800501e-15
 -3.552713678800501e-15
 -1.7763568394002505e-15
 -3.552713678800501e-15
  1.0658141036401503e-14
  1.2434497875801753e-14
  0.0

```

---

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [June 24, 2025, 4:55pm UTC](https://discourse.julialang.org/t/matrix-multiplication-and-vector-multiplication-give-different-results/130177/2 "2025-06-24T16:55:15Z")

</div>

> [@PSA: floating-point arithmetic](https://discourse.julialang.org/t/psa-floating-point-arithmetic/8678):
>
> Sometimes people are surprised by the results of floating-point calculations such as julia\> 5/6 0.8333333333333334 # shouldn't the last digit be 3? julia\> 2.6 - 0.7 - 1.9 2.220446049250313e-16 # shouldn't the answer be 0? These are not bugs in Julia. They’re consequences of the IEEE-standard 64-bit binary representation of floating-point numbers that is burned into computer hardware, which Julia and many other languages use by default. Brief explanation You can t…

---

<div class="post-metadata">

### Author: ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)
#### Post date: [June 24, 2025, 4:57pm UTC](https://discourse.julialang.org/t/matrix-multiplication-and-vector-multiplication-give-different-results/130177/3 "2025-06-24T16:57:13Z")

</div>

It’s not the problem of floating point resolution. It’s a matter of _ **different** _ implementations in matrix-matrix vs matrix-vector multiplication!

---

<div class="post-metadata">

### Author: ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [June 24, 2025, 5:03pm UTC](https://discourse.julialang.org/t/matrix-multiplication-and-vector-multiplication-give-different-results/130177/4 "2025-06-24T17:03:52Z")

</div>

Yes, this is covered in that PSA a bit further down:

> [@PSA: floating-point arithmetic](https://discourse.julialang.org/t/psa-floating-point-arithmetic/8678/8):
>
> Note there are still some ways you might get non-deterministic results:
> 
> - BLAS or LAPACK operations
> - `@fastmath`: which lets the compiler do a _lot_ of manipulations
> - `@simd`: which allows re-association of arithmetic operations to exploit SIMD operations.
> - `muladd`: allows use of either `a*b+c` or `fma(a,b,c)` depending on which is faster. This is a tricky one, as we’re increasingly making use of it (such as the recent libm work).

With lots of stuff in hardware, doing things like vectorization or reassociation (changing the order of operations) will change results.

```julia
julia> v = rand(10000);

julia> sum(v)
5004.669545612722

julia> let s = 0.0
           for x in v
               s += x
           end
           s
       end
5004.669545612701

```

Here’s an easy way to see why this happens:

```julia
julia> (1 + 1e20) - 1e20
0.0

julia> 1 + (1e20 - 1e20)
1.0

```

Matrix multiplication does all sorts of stuff that can change the order of operations, or even use reduced / improved accuracy fma operations, etc. This will naturally cause discrepancies like those you saw.

This is why people are advised to use `isapprox` (`≈`) instead of equality:

```julia
julia> m12[:, 400] == m12col400
false

julia> m12[:, 400] ≈ m12col400
true

```

These two results are equal up to the sorts of resolution guarenteed by the algorithms we use.

---

<div class="post-metadata">

### Author: ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)
#### Post date: [June 24, 2025, 5:04pm UTC](https://discourse.julialang.org/t/matrix-multiplication-and-vector-multiplication-give-different-results/130177/5 "2025-06-24T17:04:11Z")

</div>

> [@tomtom](#):
>
> It’s a matter of _ **different** _ implementations in matrix-matrix vs matrix-vector multiplication!

I don’t see any surprise. One is a [BLAS level-2](https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms#Level_2) operation (gemv), the other one is a [BLAS level-3](https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms#Level_3) operation (gemm), they normally have different implementations in BLAS libraries.
