# Different results for cholesky decomposition between 0.6.2 and 0.6.3

**URL:** <https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340>\
**Category:** General Usage\
**Tags:** linearalgebra\
**Created:** [June 1, 2018, 2:37pm UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340 "2018-06-01T14:37:35Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![matthieu](https://avatars.discourse-cdn.com/v4/letter/m/da6949/32.png) [@matthieu](https://discourse.julialang.org/u/matthieu)\
**Post date:** [June 1, 2018, 2:37pm UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340/1 "2018-06-01T14:37:35Z")

</div>

I use the following code to select a set of linearly independent columns from a matrix

```julia
function give_independent(X)
   chol = cholfact!(X'X, :U, Val{true})
   ipermute!(diag(chol.factors) .> 0, chol.piv)
end 

```

However I obtain different results in Julia 0.6.2 and 0.6.3

```julia
v1 = [1.51522, 1.91087, 1.6913, 1.88696, 1.26522, 3.19783, 3.09565, 3.03043, 4.13913, 2.27174, 1.61522, 1.43913, 2.02826, 2.45652, 0.906522, 0.543478, -0.141304, 0.528261, 5.53261, 2.98043, 2.92609, -1.87174, 0.941304, -1.53696, -1.21087, -5.63696, -8.33913, -4.82826, -8.2, -0.734783, -3.18478, -3.88913, -3.9087, -0.0130435, -1.13478, -1.10217, 2.59565, 0.530435, -0.0608696, -2.02826, -1.78478, -2.46087, 3.12826, 1.15652, -0.493478, -0.956522, 1.6587, 0.728261, 0.0326087, 3.18043, -0.973913, 0.128261, 0.741304, 2.16304, 2.28913, -6.53696, -5.03913, -13.7283, -1.2, -7.23478, -0.0847826, -0.58913, -0.908696, 0.686957, 0.0652174, -1.50217, -2.60435, 0.130435, 0.23913, 4.07174, 4.61522, 3.83913, 4.02826, 2.95652, 3.40652, 1.94348, 1.4587, 1.42826, 1.43261, 1.98043, 0.526087, 0.928261, 0.441304, -0.936957, 1.78913, -0.136957, -2.93913, -2.72826, -6.1, 4.36522]
v2 = [-25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -25.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -23.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261, -22.8261]
give_independent(hcat(v1, v1, v2))
# Julia 0.6.2. returns [false, true, true]
# Julia 0.6.3. returns [true, true , true]

```

Julia 0.6.3. does not return the result I expected since the first and second columns are linearly dependent.  
What changed between the two versions? What should I do to obtain the correct set of linearly independent vectors in 0.6.3?

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [June 1, 2018, 7:41pm UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340/2 "2018-06-01T19:41:12Z")

</div>

Thanks for the report. What’s your `versioninfo()`?

---

<div class="post-metadata">

**Author:** ![matthieu](https://avatars.discourse-cdn.com/v4/letter/m/da6949/32.png) [@matthieu](https://discourse.julialang.org/u/matthieu)\
**Post date:** [June 1, 2018, 8:23pm UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340/3 "2018-06-01T20:23:02Z")

</div>

Thanks. Here is the versioninfo()

```julia
Julia Version 0.6.2
Commit d386e40c17 (2017-12-13 18:08 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i5-7600 CPU @ 3.50GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Prescott)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, broadwell)

```

and

```julia
Julia Version 0.6.3
Commit d55cadc350 (2018-05-28 20:20 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i5-7600 CPU @ 3.50GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, broadwell)

```

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [June 1, 2018, 10:41pm UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340/4 "2018-06-01T22:41:29Z")

</div>

Interesting. This seems to be a CPU-architecture-dependent issue; I cannot reproduce on my old westmere processor. I wonder if it’s related to the different BLAS configurations. Are you using the distributed binaries or building from source?

Could you file an issue at [https://github.com/JuliaLang/julia/issues/new](https://github.com/JuliaLang/julia/issues/new) so this doesn’t get lost?

---

<div class="post-metadata">

**Author:** ![matthieu](https://avatars.discourse-cdn.com/v4/letter/m/da6949/32.png) [@matthieu](https://discourse.julialang.org/u/matthieu)\
**Post date:** [June 1, 2018, 11:57pm UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340/5 "2018-06-01T23:57:51Z")

</div>

Ok, I filed an issue there. We can continue the discussion there, but I am using distributed binaries in both cases.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [June 2, 2018, 2:38am UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340/6 "2018-06-02T02:38:07Z")

</div>

TL;DR: Small changes in the BLAS library will effect the order of operations, so there is no guarantee of identical results between minor Julia versions here.

Although your input matrix is exactly rank deficient, after a few manipulations it is perturbed, so we arrive at

> [@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…

By forming `A'A`, you make the situation worse (an almost-dependent vector could be obscured). Much more robust would be

```julia-auto
F=qrfact(A,Val{true})
ipermute!(abs.(diag(F[:R])) .> eps(10*size(A,1)*F[:R][1,1]), F[:p])

```

(The 10 is a fairly conservative tolerance factor.) Even this will fail in unusual cases, so you should check the effective rank with an SVD to be safe.

---

<div class="post-metadata">

**Author:** ![matthieu](https://avatars.discourse-cdn.com/v4/letter/m/da6949/32.png) [@matthieu](https://discourse.julialang.org/u/matthieu)\
**Post date:** [June 2, 2018, 6:17pm UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340/7 "2018-06-02T18:17:26Z")

</div>

Thanks for the great answer.  
My issue is that I do not want to use QR. The way I typically use this function is to get the set of linearly dependent vectors from a list of a few vectors, say less than 10, with potentially a huge number of elements, say 100 millions. In this case, the matrix X’X is very small. I can compute it without even forming the matrix with columns of vectors. Forming the matrix X and applying QR on it would require much more memory.  
Is there a not-so-bad way to still use cholfact? What about something that squares your condition?

```julia
    ipermute!(diag(chol.factors) .> 100 * size(X, 1)^2 * eps(chol.factors[1]), chol.piv)

```

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [June 2, 2018, 7:13pm UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340/8 "2018-06-02T19:13:54Z")

</div>

You could protect against roundoff error by using a controlled sum algorithm (similar to `sum_kbn`) to compute elements of `X'X`, then converting to `BigFloat` before factoring. (If the matrix is small, it’s less trouble to just cheat on the numerical analysis.)

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [June 3, 2018, 5:47pm UTC](https://discourse.julialang.org/t/different-results-for-cholesky-decomposition-between-0-6-2-and-0-6-3/11340/9 "2018-06-03T17:47:54Z")

</div>

10 vectors of dimension 100 million is 1e09 numbers. At double precision, with 8 bytes per number, that’s about 10e09 bytes, or 10 gigabytes. Not that big. You should be able to do a reduced QR on the full 10 x 10million data set X on a machine with 32 GB memory, or easily on a machine with 64 GB.