# Resolve sign ambiguity in SVD and Eigendecomposition

**URL:** <https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484>\
**Category:** General Usage\
**Tags:** question, linearalgebra\
**Created:** [May 19, 2021, 11:55pm UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484 "2021-05-19T23:55:44Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 19, 2021, 11:55pm UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/1 "2021-05-19T23:55:44Z")

</div>

Hello,

Given a matrix A \in \mathbb{R}^{m\times n}, my goal is to recover the SVD decomposition of A from the eigendecomposition of A^\top A and A^\top A.

A compact SVD decomposition of matrix A \in \mathbb{R}^{m\times n} gives:  
A = U\_{svd}\Sigma V\_{svd}^\top, in which \Sigma is a r \times r diagonal matrix where r = min\{m, n\}. U\_{svd} \in \mathbb{R}^{m \times r} and V\_{svd} \in \mathbb{R}^{n \times r} verify U\_{svd} U\_{svd}^\top = V\_{svd} V\_{svd}^\top = I\_r.

A is not directly accessible, I only know C\_x = A^\top A \in \mathbb{R}^{m \times m} and  
C\_y = AA^\top \in \mathbb{R}^{n \times n}.

From the spectral theorem, we know that:  
C\_x = V \Lambda\_x V^\top, C\_y = U \Lambda\_y U^\top.

If we assume that m \leq n, the eigenvalues \Lambda\_x equate the first r values of \Lambda\_y, which are the square of the singular values \Sigma.

However, we have a sign ambiguity in the columns of the different eigenvectors of C\_x and C\_y.

Therefore, U\sqrt{\Lambda\_y}V[:,1:r]^\top does not give the SVD decomposition of A! The sign of the columns of U are V are not chosen such that AV = U and A^\top U = V, as in the SVD decomposition ([Singular value decomposition - Wikipedia](https://en.wikipedia.org/wiki/Singular_value_decomposition)).

Is there a way to flip the sign of the columns of U and V to recover the SVD decomposition of A?

A simple example:

A norm of 2 for one of the columns means that they are opposite of each other.

```julia
m = 10
n = 20
r = min(m, n)
A = randn(m, n)

Cx = A'*A
Cy = A*A'

Ex = eigen(Cx; sortby = λ -> -λ)
Ey = eigen(Cy; sortby = λ -> -λ)
svdA = svd(A)

@show norm(svdA.S - sqrt.(Ex.values[1:r]))
@show norm(svdA.S - sqrt.(Ey.values[1:r]))

@show norm.(eachcol(Ex.vectors[:,1:r]-svdA.V))
@show norm.(eachcol(Ey.vectors-svdA.U));

```

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [May 20, 2021, 7:21am UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/2 "2021-05-20T07:21:08Z")

</div>

The sign information that you need is lost when you start from A^TA instead of A. (For example, you cannot distinguish between A = [1] and A = [-1] if you only know that A^TA = [1].)

---

<div class="post-metadata">

**Author:** ![fph](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fph/32/17159_2.png) [@fph](https://discourse.julialang.org/u/fph)\
**Post date:** [May 20, 2021, 7:42am UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/3 "2021-05-20T07:42:41Z")

</div>

Apart from this issue, I would advise against using this method to compute an SVD in anything serious. There are other major issues:

- worse accuracy (trouble now starts when \sigma\_1 / \sigma\_n \approx `sqrt(eps)` rather than \sigma\_1 / \sigma\_n \approx `eps`)
- even more non-uniqueness issues if there are repeated singular values (take A orthogonal, for instance).

In general you should regard computing A^TA a bit like calling `inv(A)`: 99% of the time it is a mistake, and if you shuffle your computations around there are better ways to solve the same problem.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 20, 2021, 4:00pm UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/4 "2021-05-20T16:00:16Z")

</div>

Thank you for your feedback. Annoying this loss of phase. I totally agree for the conditioning number. Let me give you a bit more context.

I would like to find the important directions in the input and output spaces for a nonlinear function f: \mathbb{R}^m \xrightarrow{} \mathbb{R}^n.

For a linear function f(x) = Ax, we can just look at the SVD decomposition of A.

To generalize this to the nonlinear case, we can find these directions by looking at the eigenvalue decomposition of the integrated inner and outer product of the Jacobian of the function f integrated over the distribution \mu:

C\_x = \int (\nabla\_x f(x))^\top (\nabla\_x f(x)) d\mu(x) \in \mathbb{R}^{m,m} and C\_y = \int (\nabla\_x f(x)) (\nabla\_x f(x))^\top d\mu(x) \in \mathbb{R}^{n,n}.

In practice, we can perform a Monte Carlo appproximation of C\_x, C\_y given M samples x^i \sim \mu:

C\_x \approx \frac{1}{M}\sum\_{i=1}^M (\nabla\_x f(x^i))^\top (\nabla\_x f(x^i)) and similarly for C\_y.

Since C\_x is positive definite, an eigendecomposition gives us:

C\_x = V \Lambda\_x V^T with V\in \mathbb{R}^{m,m}, where V an orthonormal basis for the input space. We can truncate this basis to identify the most important directions based on the decay of energy spectrum of \Lambda\_x.

Similarly, we write C\_y = U \Lambda\_y U^T with U\in \mathbb{R}^{n,n}. U is an orthonormal basis for the output space.

I don’t know how to avoid the construction and the eigendecomposition of the C\_x, C\_y matrices to find the bases U, V.

If f:\mathbb{R}^m \xrightarrow{} \mathbb{R}, then we can perform an SVD decomposition of the matrices \sqrt{C\_y} = \frac{1}{\sqrt{M}} \left[\nabla f(x^1), \ldots, \nabla f(x^M) \right] \in \mathbb{R}^{n,M} and \sqrt{C\_x} = \frac{1}{\sqrt{M}} \left[\nabla f(x^1)^\top, \ldots, \nabla f(x^M)^\top \right]\in \mathbb{R}^{m, M}, and extract the left singular vectors to get U and V.

We have C\_x = \sqrt{C\_x} \sqrt{C\_x}^\top and C\_y = \sqrt{C\_y} \sqrt{C\_y}^\top.

However, for f:\mathbb{R}^m \xrightarrow{} \mathbb{R}^n, \nabla f \in \mathbb{R}^{n,m}, and we would have to assemble matrices of (n, mM) and (m, nM), which requires too much storage.

Is there a better way to do this computation without assembling these large matrices? I can only afford finite-difference to compute the Jacobian of the function f.

---

<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:** [May 20, 2021, 5:15pm UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/5 "2021-05-20T17:15:37Z")

</div>

> [@mleprovost](#):
>
> Is there a better way to do this computation without assembling these large matrices? I can only afford finite-difference to compute the Jacobian of the function f.

Since you only want to know the “important directions”, maybe you only need the singular vectors corresponding to largest few singular values?

In that case you would be able to use iterative SVD methods (e.g. [SVDL · IterativeSolvers.jl](https://julialinearalgebra.github.io/IterativeSolvers.jl/stable/svd/svdl/) or [GitHub - Jutho/KrylovKit.jl: Krylov methods for linear problems, eigenvalues, singular values and matrix functions](https://github.com/Jutho/KrylovKit.jl) or [GitHub - JuliaLinearAlgebra/RandomizedLinAlg.jl: Randomized algorithms for numerical linear algebra in Julia](https://github.com/JuliaLinearAlgebra/RandomizedLinAlg.jl) …). Then you only need to provide a function to multiply the Jacobian by an arbitrary vector, which requires a single finite-difference operation.

---

<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:** [May 20, 2021, 5:21pm UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/6 "2021-05-20T17:21:49Z")

</div>

To answer your original question, you only need to compute the eigen-decomposition of C\_x or C\_y, not both — you can get the eigenvectors of one (for nonzero eigenvalues) from the eigenvectors of the other, and the sign/phase is then consistent for the SVD. In fact, this is one way to derive the existence of the SVD from the spectral theorem for Hermitian matrices.

For example, given a nonzero eigenvalue \lambda\_k = \sigma\_k^2 \> 0 of C\_x = A^\* A and a corresponding (orthonormal) eigenvector v\_k, you can easily show that u\_k = A v\_k / \sigma\_k is an (orthonormal) eigenvector of C\_y = A A^\* with the same eigenvalue.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 20, 2021, 7:08pm UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/7 "2021-05-20T19:08:17Z")

</div>

Thank you @stevenj for your ideas!

I tried to implement it, for the linear case f(x) = Ax, so \nabla f(x) = A, but I get some issues with the `svdl` routine from ‘IterativeSolvers.jl’

```julia
using LinearAlgebra
using LinearMaps
using IterativeSolvers

Nx = 20
Ny = 10

A = randn(Ny, Nx)
Ne = 500
X = randn(Nx, Ne)

function mapsqrtCx!(Ny, Ne, vout, vin)
    fill!(vout, 0.0)
    
    @inbounds for j=1:Ne
        # chunk of dimension Ny
        vj = view(vin, (j-1)*Ny+1:j*Ny)
        # A will be the Jacobian of the function f evaluated at x^i
        mul!(vout, A', vj, 1.0, 1.0)
    end
    vout .*= 1/sqrt(Ne)
end

function mapsqrtCxtranspose!(Ny, Ne, vout, vin)
    fill!(vout, 0.0)

    @inbounds for j=1:Ne
        # chunk of dimension Ny
        vj = view(vout, (j-1)*Ny+1:j*Ny)
        # A will be the Jacobian of the function f evaluated at x^i
        mul!(vj, A, vin, 1.0, 1.0)
    end
    vout .*= 1/sqrt(Ne)
end

function mapsqrtCy!(Nx, Ne, vout, vin)
    fill!(vout, 0.0)
    
    @inbounds for j=1:Ne
        # chunk of dimension Nx
        vj = view(vin, (j-1)*Nx+1:j*Nx)
        # A will be the Jacobian of the function f evaluated at x^i
        mul!(vout, A, vj, 1.0, 1.0)
    end
    vout .*= 1/sqrt(Ne)
end

function mapsqrtCytranspose!(Nx, Ne, vout, vin)
    fill!(vout, 0.0)

    @inbounds for j=1:Ne
        # chunk of dimension Nx
        vj = view(vout, (j-1)*Nx+1:j*Nx)
        # A will be the Jacobian of the function f evaluated at x^i
        mul!(vj, A', vin, 1.0, 1.0)
    end
    vout .*= 1/sqrt(Ne)
end

sqrtCx = LinearMap((vout, vin) -> mapsqrtCx!(Ny, Ne, vout, vin), 
                   (vout, vin) -> mapsqrtCxtranspose!(Ny, Ne, vout, vin), 
                   Nx, Ny*Ne; ismutating = true)

sqrtCy = LinearMap((vout, vin) -> mapsqrtCy!(Nx, Ne, vout, vin), 
                   (vout, vin) -> mapsqrtCytranspose!(Nx, Ne, vout, vin), 
                   Ny, Nx*Ne; ismutating = true)

cache = ones(Ny, Ne)
@time 1/sqrt(Ne)*sum(A'*cache;dims = 2)

v = ones(Ny*Ne);
@time sqrtCx*v

@show norm(1/sqrt(Ne)*sum(A'*ones(Ny, Ne);dims = 2)-sqrtCx*v)

```

12.892 μs (10 allocations: 39.56 KiB)  
40.500 μs (501 allocations: 23.59 KiB)  
norm((1 / sqrt(Ne)) \* sum(A \* ones(Nx, Ne); dims = 2) - sqrtCx \* v) = 7.149791431641379e-12

Any ideas to reduce the computational time for `mapsqrtCx!`?

Now, we can use the routine `svdl` from IterativeSolvers.jI:

```julia
maxrank = 5
@time SCx, errorCx = svdl(sqrtCx; nsv = maxrank, vecs = :left);
@time SCy, errorCy = svdl(sqrtCy; nsv = maxrank, vecs = :left);

svdA = svd(A)
# check singular values
@show norm(SCx.S[1:maxrank] - svdA.S[1:maxrank])
@show norm(SCy.S[1:maxrank] - svdA.S[1:maxrank])

```

---

<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:** [May 20, 2021, 9:09pm UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/8 "2021-05-20T21:09:01Z")

</div>

> [@mleprovost](#):
>
> Any ideas to reduce the computational time for `g!` ?

matrix–matrix multiplications will almost always be much more efficient than performing the same operation as a loop of matrix–vector multiplications. A factor of 4 isn’t bad.

In any case, before you wild on micro-optimizations, in your real application won’t the dominant cost be in evaluating your nonlinear function f(x)?

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 20, 2021, 9:21pm UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/9 "2021-05-20T21:21:36Z")

</div>

Yes, the dominant cost is clearly to evaluate/differentiate f. I think the problem is that we need multiple evaluation of `sqrtCx v` for this iterative solver, therefore we need to evaluate the Jacobian of the nonlinear function f for all the different samples for a single evaluation of the product `sqrtCx v`.

I forgot the transpose version of `g` to compute the left singular vectors of `sqrtCx`.

I am not very familiar with iterative SVD solvers, but I think this is working correctly as long as the spectral gap between two consecutive singular values is big enough. On this toy problem, I have found that the singular vectors are correct if I stop at the 4th value, but stopping at the 5th value gives a very poor prediction of the singular left vectors…

---

<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:** [May 21, 2021, 1:17am UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/10 "2021-05-21T01:17:17Z")

</div>

> [@mleprovost](#):
>
> therefore we need to evaluate the Jacobian of the nonlinear function f ff for all the different samples for a single evaluation of the product `sqrtCx v`

You are not evaluating the whole Jacobian matrix, right? For a Jacobian–vector product you only need a directional derivative, e.g. by a single finite difference.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 21, 2021, 4:47am UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/11 "2021-05-21T04:47:55Z")

</div>

I am probably missing something important here.

The input vector has a dimension m \sim 100, and the output has dimension n \sim 20-40, and I have about M \sim 100-500 samples.  
So far, I was computing the entire Jacobian matrix for the different ensemble members and storing it.  
Do you mean that we can compute the gradient of f(x)^\top a with respect to x for \sqrt{C\_x} to get \nabla f(x)^\top a, where a \in \mathbb{R}^n?

For f:\mathbb{R}^m \xrightarrow{} \mathbb{R}^n, let us define the matrices \sqrt{C\_x} = \frac{1}{\sqrt{M}} \left[\nabla f(x^1)^\top, \ldots, \nabla f(x^M)^\top \right]\in \mathbb{R}^{m, nM} and \sqrt{C\_y} = \frac{1}{\sqrt{M}} \left[\nabla f(x^1), \ldots, \nabla f(x^M) \right] \in \mathbb{R}^{n,mM}.

We have C\_x = \sqrt{C\_x} \sqrt{C\_x}^\top and C\_y = \sqrt{C\_y} \sqrt{C\_y}^\top.

I don’t know how to do to get \nabla f(x)b, (with b \in \mathbb{R}^m) as the gradient of something with respect to x?  
Do you have some references on these techniques?

For these iterative solvers, do we know how many iterations they need to converge for the leading singular values?

In the linear case, based on the characterization of the SVD that you wrote above, we have:

\sqrt{C\_x} \begin{bmatrix} u\_i\\ \vdots \\ u\_i \end{bmatrix} = \sqrt{M} \sqrt{\sigma\_i} v\_i.

I believe this doesn’t hold when f is a nonlinear function? Meaning that we need to perform two SVD decompositions : one for \sqrt{C\_x} and one for \sqrt{C\_y}.

---

<div class="post-metadata">

**Author:** ![ettersi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ettersi/32/6829_2.png) [@ettersi](https://discourse.julialang.org/u/ettersi)\
**Post date:** [May 21, 2021, 7:46am UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/13 "2021-05-21T07:46:56Z")

</div>

> [@mleprovost](#):
>
> C\_x = \int (\nabla\_x f(x))^\top (\nabla\_x f(x)) d\mu(x) \in \mathbb{R}^{m,m} and C\_y = \int (\nabla\_x f(x)) (\nabla\_x f(x))^\top d\mu(x) \in \mathbb{R}^{n,n}.

I’m not sure why you do (\nabla\_x f(x))^\top (\nabla\_x f(x)) and (\nabla\_x f(x)) (\nabla\_x f(x))^\top in these integrals. Wouldn’t an SVD of

\int \nabla f(x) \, d\mu(x)

give you similar information without the problem of squared condition numbers? (I’m no expert on this part, so this is a genuine question.)

> [@mleprovost](#):
>
> I don’t know how to do to get \nabla f(x)b, (with b \in \mathbb{R}^m) as the gradient of something with respect to x? Do you have some references on these techniques?

If the gradient exists, then by definition we have

\nabla f(x) \, b = \lim\_{\varepsilon \to 0} \frac{f(x + \varepsilon b) - f(x)}{\varepsilon}.

Thus, you can compute \nabla f(x) \, b with a single application of finite differences to \varepsilon \mapsto f(x + \varepsilon b).

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [May 21, 2021, 5:35pm UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/15 "2021-05-21T17:35:59Z")

</div>

That’s a good question, a proper justification of these Gramian matrices is described in this paper:

> **[Data-Free Likelihood-Informed Dimension Reduction of Bayesian Inverse Problems](https://arxiv.org/abs/2102.13245)**
>
> Identifying a low-dimensional informed parameter subspace offers a viable path to alleviating the dimensionality challenge in the sampled-based solution to large-scale Bayesian inverse problems. This paper introduces a novel gradient-based dimension...

You can also look for active subspaces  
[https://epubs.siam.org/doi/book/10.1137/1.9781611973860?mobileUi=0&](https://epubs.siam.org/doi/book/10.1137/1.9781611973860?mobileUi=0&)

For example in 2D, if f(x\_1, x\_2) = [x\_1; x\_2^2], so \partial f\_2/\partial x\_2 = 2 x\_2. If you take the expectation of the Jacobian over a zero-mean symmetric distribution along the x\_2 axis (e.g. a standard Gaussian distribution), then the expectation of this component is clearly zero, and the most informative direction in the input space is missed.

Thank you for your explanation, how would you compute \nabla f(x)^\top a with a single application of finite differences? We can say that this is the gradient of f(x)^\top a with respect to x, but I don’t think that we can use the trick of the directional derivative as for \nabla f(x) b?

---

<div class="post-metadata">

**Author:** ![ettersi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ettersi/32/6829_2.png) [@ettersi](https://discourse.julialang.org/u/ettersi)\
**Post date:** [May 22, 2021, 11:46am UTC](https://discourse.julialang.org/t/resolve-sign-ambiguity-in-svd-and-eigendecomposition/61484/16 "2021-05-22T11:46:29Z")

</div>

> [@mleprovost](#):
>
> For example in 2D, if f(x\_1, x\_2) = [x\_1; x\_2^2] f(x1,x2)=[x1;x22]f(x\_1, x\_2) = [x\_1; x\_2^2] , so \partial f\_2/\partial x\_2 = 2 x\_2 ∂f2/∂x2=2x2\partial f\_2/\partial x\_2 = 2 x\_2 . If you take the expectation of the Jacobian over a zero-mean symmetric distribution along the x\_2 x2x\_2 axis (e.g. a standard Gaussian distribution), then the expectation of this component is clearly zero, and the most informative direction in the input space is missed.

Makes sense.

> [@mleprovost](#):
>
> how would you compute \nabla f(x)^\top a ∇f(x)⊤a\nabla f(x)^\top a with a single application of finite differences? We can say that this is the gradient of f(x)^\top a f(x)⊤af(x)^\top a with respect to x xx , but I don’t think that we can use the trick of the directional derivative as for \nabla f(x) b ∇f(x)b\nabla f(x) b ?

I believe the only way to evaluate gradients through finite differencing is to compute the derivatives with respect to each input dimension. If you are open to replacing finite differences with automatic differentiation, then what you would have to do is to replace forward mode AD with reverse mode AD.