# Floating number matrix, equal matrix in whole, but not equal column

**URL:** <https://discourse.julialang.org/t/floating-number-matrix-equal-matrix-in-whole-but-not-equal-column/101200>\
**Category:** Numerics\
**Created:** [July 5, 2023, 9:39am UTC](https://discourse.julialang.org/t/floating-number-matrix-equal-matrix-in-whole-but-not-equal-column/101200 "2023-07-05T09:39:37Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![maxchendt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxchendt/32/42976_2.png) [@maxchendt](https://discourse.julialang.org/u/maxchendt)\
**Post date:** [July 5, 2023, 9:39am UTC](https://discourse.julialang.org/t/floating-number-matrix-equal-matrix-in-whole-but-not-equal-column/101200/1 "2023-07-05T09:39:37Z")

</div>

`norm(Y - B * A, Inf)` is 0, but when comparing individual columns, such as `Y[:,4] - B * A[:,4]` is not zero. Why?

```julia
using TranscodingStreams, CodecXz, Serialization, Downloads
using LinearAlgebra

function urlData(url, fn)
    tFile = joinpath(tempdir(), fn)
    if !isfile(tFile)
        Downloads.download(url, tFile)
    end
    #xz = last(fn, 3)
    io = open(tFile, "r")
    if lowercase(last(fn, 3)) == ".xz"
        io = TranscodingStream(XzDecompressor(), io)
    end
    t = deserialize(io)
    close(io)
    return t
end

function main()

    url = "https://github.com/maxchendt/FuturesHedge/raw/master/data/YA4.jls.xz"
    fn = "YA4.jls.xz"
    B, Y, A = urlData(url, fn)

    #display(Y - B * A)
    display(norm(Y - B * A, Inf))
    display(Y[:,4] - B * A[:,4])
    #y = B * A[:,4]
    #display(Y[:,4] - y)

    nothing

end

main()
nothing

```

the output

```julia
0.0
3-element Vector{Float64}:
 -1.3877787807814457e-17
  0.0
  2.7755575615628914e-17

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [July 5, 2023, 10:59am UTC](https://discourse.julialang.org/t/floating-number-matrix-equal-matrix-in-whole-but-not-equal-column/101200/2 "2023-07-05T10:59:23Z")

</div>

That is probably due to `norm` being optimized and performing the floating point operations in a slightly different order than you would do column-by-column. If you test `isapprox(Y[:,4], B * A[:,4])` to disregard this kind of numerical errors, you will get `true`

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [July 5, 2023, 11:53am UTC](https://discourse.julialang.org/t/floating-number-matrix-equal-matrix-in-whole-but-not-equal-column/101200/3 "2023-07-05T11:53:41Z")

</div>

> [@maxchendt](#):
>
> `norm(Y - B * A, Inf)` is 0, but when comparing individual columns, such as `Y[:,4] - B * A[:,4]` is not zero. Why?

Due to rounding errors, `B * A` gives only _approximately_ the same result as multiplying `B` separately by each column of `A`.

```julia
julia> A, B = randn(1000,1000), randn(1000,1000);

julia> (B * A)[:,4] == B * A[:,4]
false

julia> (B * A)[:,4] ≈ B * A[:,4]
true

```

The reason is that `B * A` uses a different (faster) algorithm for multiplying the matrices than column-by-column, which does the operations in a different order, and floating-point arithmetic is not associative. Even just summing the same numbers in a different order produces a different result due to roundoff errors accumulating differently.

---

<div class="post-metadata">

**Author:** ![maxchendt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maxchendt/32/42976_2.png) [@maxchendt](https://discourse.julialang.org/u/maxchendt)\
**Post date:** [July 5, 2023, 10:14pm UTC](https://discourse.julialang.org/t/floating-number-matrix-equal-matrix-in-whole-but-not-equal-column/101200/4 "2023-07-05T22:14:55Z")

</div>

> [@stevengj](#):
>
> floating-point arithmetic is not associative

Thanks so much!

I have another big matrix, the difference is 10^-8. I double check and check my codes, now I am not puzzled.