# SVD gives different results from Python with MWE

**URL:** <https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674>\
**Category:** General Usage\
**Tags:** linearalgebra, python\
**Created:** [March 7, 2021, 1:17pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674 "2021-03-07T13:17:03Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![kadir-gunel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kadir-gunel/32/35790_2.png) [@kadir-gunel](https://discourse.julialang.org/u/kadir-gunel)\
**Post date:** [March 7, 2021, 1:17pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/1 "2021-03-07T13:17:03Z")

</div>

Hello,

I have a python written code, and it uses svd from numpy. And I am trying to port the same code in Julia.  
The problem I am facing is that after some point, even though I read the same files for both programming languages, svd method in two languages decomposes the same matrix differently. Hence, solution changes and I cannot replicate the experiment.

I have multiple svd usages inside the code, the first two gives nearly the exactly decomposition. (Everything is in Float32) But after the third svd function call my results and the original(python) code starts to differ. The only difference from the first 2 svd functions is that the third one removes the ` full_matrices=False` flag in Python code.

```python
 u, s, vt = np.linalg.svd(xsim[:idx], full_matrices=False) # get the same result
 u, s, vt = np.linalg.svd(zsim[:idx], full_matrices=False) # get the same result
 u, s, vt = np.linalg.svd(z[trg_indices].T.dot(x[src_indices])) # third call to svd; where things get messy.

```

On the other hand, in Julia, even thoug I changed the algorithm to QRIteration , and full to true. I am not getting the same results.

I did some research and both numpy and Julia uses LAPACK package underneath. Then what should be the correct way to get the same results ?

This is what I get when I write `print(numpy.show_config())`

```python
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]

```

For julia could someone please help me if all set is the same as Python ?

B.R.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [March 7, 2021, 1:25pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/2 "2021-03-07T13:25:26Z")

</div>

It would be much easier to help you if you could create a MWE that is give us the specific matrices `xsim[:idx]`, `zsim[:idx]`, and `z[trg_indices].T.dot(x[src_indices])` such that we can reproduce the behaviour.

On a general note: Have you taken account for the fact that Julia uses column-major ordering and Python row-major ordering? That is matrices in Julia in Python might differ by a transpose operation.

---

<div class="post-metadata">

**Author:** ![kadir-gunel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kadir-gunel/32/35790_2.png) [@kadir-gunel](https://discourse.julialang.org/u/kadir-gunel)\
**Post date:** [March 7, 2021, 1:27pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/3 "2021-03-07T13:27:35Z")

</div>

Yeah ,you are right. But those matrices are huge files nearly (one is 1.1GB and the other is 512MB). If you want I can share you the links.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [March 7, 2021, 1:29pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/4 "2021-03-07T13:29:44Z")

</div>

Ideally, you should try to create a **minimal** example. For example, If you only consider the submatrix, say `M[1:100,1:100]`, do Numpy and Julia still disagree?

Having said that, sharing a link to a large file is still better than sharing nothing 🙂

---

<div class="post-metadata">

**Author:** ![kadir-gunel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kadir-gunel/32/35790_2.png) [@kadir-gunel](https://discourse.julialang.org/u/kadir-gunel)\
**Post date:** [March 7, 2021, 1:33pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/6 "2021-03-07T13:33:14Z")

</div>

> [@carstenbauer](#):
>
> Ideally, you should try to create a **minimal** example. For example, If you only consider the submatrix, say `M[1:100,1:100]` , do Numpy and Julia still disagree?

Actually I am reading submatrices of these huge matrices. The first 4000 words.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [March 7, 2021, 1:36pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/7 "2021-03-07T13:36:22Z")

</div>

Instead of posting a link to a “raw” file that then needs to get loaded and processed by a function you could make it easier for the other side and dump the final matrix, i.e. `xsim[:idx]` to a reasonable file format that can be read directly. For example, using [NPZ.jl](https://github.com/fhs/NPZ.jl) which is compatible with `numpy.save` and `numpy.load`.

This way you can simply load the file in Julia and Python (i.e. numpy) and assure that you actually operate on precisely the same data.

(You want to abstract away your context as much as possible when creating a MWE.)

---

<div class="post-metadata">

**Author:** ![kadir-gunel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kadir-gunel/32/35790_2.png) [@kadir-gunel](https://discourse.julialang.org/u/kadir-gunel)\
**Post date:** [March 7, 2021, 1:38pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/8 "2021-03-07T13:38:50Z")

</div>

> [@carstenbauer](#):
>
> Instead of posting a link to a “raw” file that then needs to get loaded and processed by a function you could make it easier for the other side and dump the final matrix, i.e. `xsim[:idx]` to a reasonable file format that can be read directly. For example, using [NPZ.jl](https://github.com/fhs/NPZ.jl) which is compatible with `numpy.save` and `numpy.load` .

Thank didn’t know that. I will upload them.

---

<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:** [March 7, 2021, 1:39pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/9 "2021-03-07T13:39:49Z")

</div>

> [@kadir-gunel](#):
>
> I have multiple svd usages inside the code, the first two gives nearly the exactly decomposition. (Everything is in Float32) But after the third svd function call my results and the original(python) code starts to differ.

Differ in what way? Realize that you will often get slightly different results due to roundoff errors, and there is also some arbitrariness in the choice of phase of the singular vectors (or their sign, for real matrices) that can change between versions of LAPACK/BLAS. If you are using the full (not thin) SVD, then there can also be arbitrary rotations in the nullspace portions of the basis (as well as rotations in vectors for pairs of equal singular values).

---

<div class="post-metadata">

**Author:** ![kadir-gunel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kadir-gunel/32/35790_2.png) [@kadir-gunel](https://discourse.julialang.org/u/kadir-gunel)\
**Post date:** [March 7, 2021, 2:06pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/10 "2021-03-07T14:06:16Z")

</div>

Could you please try svd with these [both source and target embeddings](https://www.dropbox.com/sh/fli7mv9tqjbg5fw/AAC0e1Amlc32e93fNziyDLXVa?dl=0) ?

For the MWE:

```julia

using LinearAlgebra
using NPZ # for loading uploaded matrices

X = npzread("./path_to_file/src_emb.npy")
Z = npzread("./path_to_file/trg_emb.npy")

F = svd(Z' * X, full=true)
W = F.V * F.U'

```

W should show :

```julia
300×300 Array{Float32,2}:
 -0.136035 -0.0777119 0.0617536 0.0322536 -0.0803463 -0.0503565 … 0.095194 -0.0131295 -0.0396312 -0.00628574 0.0511491
 -0.00174201 0.0281835 0.0454273 0.0358372 0.0529457 -0.0639435 -0.0942243 0.0325997 0.000161858 -0.0514902 -0.0498046
  0.0676922 0.237949 0.0913208 -0.0106957 -0.00361676 -0.0334244 -0.0415728 0.0830503 0.0839673 -0.0716634 -0.0183218
  0.180583 -0.0958469 -0.0631438 -0.000189014 0.0190076 -0.0152651 0.0361847 -0.00151322 0.0589926 0.0360943 -0.120823
  0.0309126 0.0636542 0.0391062 -0.0415557 -0.107748 0.0503444 -0.00643703 -0.0312398 0.0658869 0.00057633 -0.0033062
  ⋮ ⋮ ⋱ ⋮                                                   
 -0.00427003 -0.0422672 0.0433433 0.0792277 -0.0729017 -0.0360674 0.0137698 0.00893876 0.0135633 -0.00463451 -0.0784995
  0.0467229 0.0159696 -0.0542639 -0.0673708 0.0118304 0.0391193 0.121941 -0.041253 -0.0230913 -0.0481352 -0.0489984
  0.0685667 0.00899705 -0.0338196 -0.102029 -0.072511 0.074575 0.0204513 0.0675888 0.0863168 -0.00807602 0.00272199
  0.0906389 0.08509 -0.0440936 0.0993734 0.0830268 0.0200159 0.0228065 0.067321 0.0181914 -0.0103052 -0.0205173

```

Whereas python:

```python

import numpy as np 

X = np.load('./path_to_/src_emb.npy')
Z = np.load('./path_to_/trg_emb.npy')

u,s,vt = np.linalg.svd(Z.T.dot(X))

w = vt.T.dot(u.T)

```

w should show

```julia
array([[ 0.02077192, -0.09222455, 0.16809827, ..., -0.07213789,
         0.01867277, -0.03105616],
       [-0.0070491 , -0.04300215, 0.02156948, ..., 0.08834429,
        -0.01833235, -0.01745731],
       [ 0.08617274, 0.15471898, 0.01507143, ..., 0.01034247,
        -0.11459561, -0.01204915],
       ...,
       [ 0.00199408, 0.01696967, -0.04052072, ..., 0.00589352,
         0.01136548, -0.05689107],
       [ 0.01693112, -0.05401053, -0.03806448, ..., 0.01929145,
        -0.03649945, 0.03389705],
       [ 0.12757824, 0.05976016, -0.03483151, ..., 0.04901506,
        -0.09714089, -0.04803994]], dtype=float32)

```

> [@stevengj](#):
>
> Realize that you will often get slightly different results due to roundoff errors,

Yes, I understand and expect that. But after decomposing the mutliplication of both matrices, all matrices U, V start similar but they are quite dissimilar at later column space. I am not claiming anything and do not know much about how things work under the hood. I am just wondering how to make both programs give the same (similar) results.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [March 7, 2021, 3:13pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/11 "2021-03-07T15:13:24Z")

</div>

Well, the results are certainly similar:

```julia
using NPZ, PyCall, Statistics, LinearAlgebra
np = pyimport("numpy");

adiff(A,B) = maximum(abs.(abs.(A) .- abs.(B)))

mat_src = npzread("src_emb.npy"); # Float32 data
U, S, Vt = np.linalg.svd(mat_src, full_matrices=false)
F = svd(mat_src);
display((adiff(F.U, U), adiff(F.S, S), adiff(F.Vt, Vt)))
display((abs.(F.U) ≈ abs.(U), F.S ≈ S, abs.(F.Vt) ≈ abs.(Vt)))

mat_trg = npzread("trg_emb.npy"); # Float32 data
U, S, Vt = np.linalg.svd(mat_trg, full_matrices=false)
F = svd(mat_trg);
display((adiff(F.U, U), adiff(F.S, S), adiff(F.Vt, Vt)))
display((abs.(F.U) ≈ abs.(U), F.S ≈ S, abs.(F.Vt) ≈ abs.(Vt)))

```

I get:

```julia
(0.0002946444f0, 9.1552734f-5, 0.0007101372f0)
(false, true, false)

(4.1430816f-5, 0.00021362305f0, 0.00010115653f0)
(true, true, true)

```

---

<div class="post-metadata">

**Author:** ![kadir-gunel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kadir-gunel/32/35790_2.png) [@kadir-gunel](https://discourse.julialang.org/u/kadir-gunel)\
**Post date:** [March 7, 2021, 3:40pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/12 "2021-03-07T15:40:28Z")

</div>

Thank you for your response.

In your sample, you decomposed individual matrices. And I assume you applied svd on individual matrices in order to show that the difference is so small that it does not matter if there is a multiplication operation (or another op.) between the matrices.  
And as you are trying to explain the differences between decomposed values are really small which are negligible.

But when I do for the multiplication of two matrices the result is :

```julia
julia> display((adiff(F.U, U), adiff(F.S, S), adiff(F.Vt, Vt)))
(0.47538608f0, 0.0078125f0, 0.6747172f0)

julia> display((abs.(F.U) ≈ abs.(U), F.S ≈ S, abs.(F.Vt) ≈ abs.(Vt)))
(false, true, false)

```

For the singular values (S and F.S), the difference again is very small. But cannot say the same for the orthogonal matrices U, V .

Could you please clarify what should I do ?  
Thank you in advance

---

<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:** [March 7, 2021, 4:51pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/13 "2021-03-07T16:51:52Z")

</div>

> [@kadir-gunel](#):
>
> `W = F.V * F.U'`

I think the basic problem here is that this `W` is an ill-conditioned function of `A`, so it will be extremely sensitive to small differences in roundoff errors (in any language).

It correspond to taking A^T = V \Sigma U^T and then forming a new matrix W = A^T U \Sigma^{-1} U^T = VU^T. This is a problem because whenever A has small singular values the matrix U \Sigma^{-1} U^T that you are (implicitly) multiplying by is ill-conditioned. (In your example, the 300\times 300 matrix A has numerical rank \approx 133 — its near-zero singular values are basically numerical noise that you are amplifying to unity).

A symptom of this is that, even in Python, you get a completely different W matrix if you do the computation in double vs. single precision:

```julia
import numpy as np 

X = np.load('src_emb.npy')
Z = np.load('trg_emb.npy')

u,s,vt = np.linalg.svd(Z.T.dot(X))
w = vt.T.dot(u.T)

u64,s64,vt64 = np.linalg.svd(np.float64(Z).T.dot(np.float64(X)))
w64 = vt64.T.dot(u64.T)

np.linalg.norm(w - w64, 2)

```

gives `1.99999...` (i.e. roughly 2, which is what you would expect for the operator norm of the difference between two _random_ unitary matrices).

You may need to re-think the computation you are performing if it depends on this matrix. (In some cases the final result may be okay if the algorithm is sufficiently clever/careful, even if intermediate matrices like W vary wildly.)

---

<div class="post-metadata">

**Author:** ![kadir-gunel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kadir-gunel/32/35790_2.png) [@kadir-gunel](https://discourse.julialang.org/u/kadir-gunel)\
**Post date:** [March 7, 2021, 5:47pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/15 "2021-03-07T17:47:55Z")

</div>

> You may need to re-think the computation you are performing if it depends on this matrix. (In some cases the final result may be okay if the algorithm is sufficiently clever/careful, even if intermediate matrices like W vary wildly.)

svd occurs inside in a loop and since W is ill-conditioned (didn’t know that) 2 programs start to differ immediately. The one in python optimizes the result much better where as the Juila code stucks and finishes early.

Thank you very much for your time and clean explanation.

Cheers from Istanbul

---

<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:** [March 7, 2021, 9:11pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/16 "2021-03-07T21:11:36Z")

</div>

PS. You can also pass `alg=LinearAlgebra.QRIteration()` to Julia’s `svd` function to use a slower (but slightly more accurate) algorithm, in case that helps. I would still be worried about relying on W=VU^\* if you want a robust algorithm, however.

---

<div class="post-metadata">

**Author:** ![kadir-gunel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kadir-gunel/32/35790_2.png) [@kadir-gunel](https://discourse.julialang.org/u/kadir-gunel)\
**Post date:** [March 8, 2021, 6:00am UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/17 "2021-03-08T06:00:44Z")

</div>

> [@stevengj](#):
>
> PS. You can also pass `alg=LinearAlgebra.QRIteration()` to Julia’s `svd` function to use a slower (but slightly more accurate) algorithm, in case that helps. I would still be worried about relying on W=VU∗W=VU^\* if you want a robust algorithm, however.

I have tried `QRIteration()` but nothing changed.

Thank you

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [March 9, 2021, 2:32pm UTC](https://discourse.julialang.org/t/svd-gives-different-results-from-python-with-mwe/56674/18 "2021-03-09T14:32:53Z")

</div>

Having exhausted other options, you can try to reverse engineer Numpy to see if there are preprocessing the matrix somehow:

> <https://github.com/numpy/numpy/blob/v1.20.0/numpy/linalg/linalg.py#L1482-L1674>

It looks like they use a subset of LAPACK processed through `f2c` in the lapack\_lite package.

> **[numpy/numpy/linalg/lapack\_lite at v1.20.0 · numpy/numpy](https://github.com/numpy/numpy/tree/v1.20.0/numpy/linalg/lapack_lite)**
>
> v1.20.0/numpy/linalg/lapack\_lite