Julia: Different Results for Cholfact decomposition between 0.6.2 and 0.6.3.

[matthieugomez]

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

function base_cols(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

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]
base_cols(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.

Here is the versioninfo()

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

This issue also causes my package to fail on Travis (while it did not fail with 0.6.2.)

[chriselrod]

I can reproduce the 0.6.3 side of things on a computer with MKL:

julia> versioninfo()
Julia Version 0.6.3
Commit d55cadc350 (2018-05-28 20:20 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i3-4010U CPU @ 1.70GHz
  WORD_SIZE: 64
  BLAS: libmkl_rt
  LAPACK: libmkl_rt
  LIBM: libimf
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

julia> base_cols(hcat(v1, v1, v2))
3-element BitArray{1}:
 true
 true
 true

I get the same answer on a 7 day old master (also with MKL) after replace ipermute! with invpermute!.

I don't have 0.6.2 anymore, but could reinstall it.

This implies the problem isn't with OpenBLAS.

[andreasnoack]

This method isn't a reliable way of detecting linear dependency. The Cholesky factorization will introduce rounding errors and the tiniest perturbation will cause the diagonal element to be non-zero. We upgraded the version of OpenBLAS from 0.2.20 to 0.3.0 between Julia 0.6.2 and 0.6.3 which probably explains the difference. However, it was probably pure chance that you got an exact zero in 0.6.2.

[atbug]

I believe this is not an issue introduced in julia 0.6.3. The following result is from JuliaPro 0.6.2.2

julia> versioninfo()
Julia Version 0.6.2
Commit d386e40* (2017-12-13 18:08 UTC)
Platform Info:
  OS: Linux (x86_64-unknown-linux-gnu)
  CPU: Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.40GHz
  WORD_SIZE: 64
  BLAS: libmkl_rt
  LAPACK: libmkl_rt
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

julia> base_cols(hcat(v1, v1, v2))
3-element BitArray{1}:
 true
 true
 true

[matthieugomez]

Thanks. What is the recommended way? Should I close this issue?

[andreasnoack]

You'd have to use a tolerance. The pivoted Cholesky allows for a tol argument which uses a default value if it is negative so you could do

F = cholfact(Symmetric(X'X), Val{true}, tol = -1.0)

and then use rank(F) to determine how many independent columns there are. Actually, I think that we should use the default value instead on 0.0 which we currently do since the primary use of the pivoted Cholesky probably is to determine the rank.

[matthieugomez]

Thanks a lot!