Singular value decomposition of complex valued matrices does not match mathematical definition

[IvanYashchuk]

📚 Documentation

Now when torch.svd supports complex input both for CPU and CUDA. Let's discuss the docs.
Currently, SVD documentation is fully in agreement with the implementation. This function returns(U, S, V) such that input = U @ diag(S) @ V.T
However, the mathematical definition of SVD is the decomposition of a matrix M such that , where V^* is the conjugate transpose of V.

What is the opinion on this? Should the documentation and implementation match math and be updated such that A = U @ diag(S) @ V.T.conj()?

Link to torch.svd documentation.

cc @ezyang @gchanan @zou3519 @bdhirsh @heitorschueroff @anjali411 @dylanbespalko @mruberry @jlin27 @vishwakftw @jianyuh @nikitaved @pearu @vincentqb

[pearu]

xref: #45063 - conjugate transpose operator feature request

[pearu]

IIUC, the current SVD documentation uses transpose because previously only real input was supported (and hence using conjugate transpose would have been unnecessary). With the new complex input support in SVD, I would suggest updating the SVD documentation to use conjugate transpose to match the mathematical definition of general SVD. While doing so, also mention that when M is real then V will be real and conjugate transpose falls back to transpose.

[ezyang]

Yes, please just generalize the docs, in this case it seems to be straightforward and the right thing to do.

[nikitaved]

The current SVD implementation returns V that is silently conjugated. I guess it makes sense to fix that too, is that correct?

[mruberry]

The current SVD implementation returns V that is silently conjugated. I guess it makes sense to fix that too, is that correct?

As far as I understand what you're saying: yes.

[antocuni]

I'd like to add my 2 cents to this discussion. By looking at the source code, it seems to me that the current behavior or torch.svd was not designed, but simply a byproduct of a bug. In particular, look at _create_U_S_VT:

pytorch/aten/src/ATen/native/LinearAlgebraUtils.h

Line 221 in a814231

static inline std::tuple<Tensor, Tensor, Tensor> _create_U_S_VT(const Tensor& input, bool some, bool compute_uv) {

As the name suggests, it should create a VT tensor, which is done here:

pytorch/aten/src/ATen/native/LinearAlgebraUtils.h

Lines 246 to 247 in a814231

	// VT should be a row-major or a batch of row-major matrices
	Tensor VT_empty;

However, the comment is wrong: the lapack's function _gesdd returns a VT in column-major order.

The rest of the code keeps calling it VT; e.g., this is _svd_helper_cpu returning a tuple (U, S, VT):

pytorch/aten/src/ATen/native/BatchLinearAlgebra.cpp

Line 1058 in a814231

return std::make_tuple(U_working_copy, S_working_copy, VT_working_copy);

However, since we are creating the tensor with the wrong strides, we are effectively returning VT.transpose(), which is V. To summarize:

• numpy.linalg.svd returns VT
• lapack retusn VT
• our implementation returns something which is called VT (but it's not)

So, it looks to me that the intention was to return VT as well, but because of the bug in _create_U_S_VT, we started to return V and instead of fixing the bug we simply "fixed" the documentation.

There isn't much that we can to at this point without breaking compatibility, but this confusing/unexpected/nonstandard behavior of svd might indicate that we should deprecate it in favor or torch.linalg.svd, once it lands.

Note: in my WIP PR #45562 I fixed _create_U_S_VT by creating VT with the correct strides, and transpose_()ing it upon return to keep b/c.
In the same PR I am also going to fix the docs for torch.svd.

[mruberry]

Sounds great, @antocuni! I agree we should consider deprecating torch.svd once torch.linalg.svd lands. It will also be confusing to have a method torch.Tensor.svd when both functions are in because users may not know which svd it refers to.

[antocuni]

It will also be confusing to have a method torch.Tensor.svd when both functions are in because users may not know which svd it refers to.

This is something which I didn't think about. What should Tensor.svd do if/when we remove torch.svd? Ideally it would do the same as torch.linalg.svd, but this will break compatibility.

In the same PR I am also going to fix the docs for torch.svd.

FWIW, this is the commit which aims to fix the documentation:
97662da

[mruberry]

It will also be confusing to have a method torch.Tensor.svd when both functions are in because users may not know which svd it refers to.

This is something which I didn't think about. What should Tensor.svd do if/when we remove torch.svd? Ideally it would do the same as torch.linalg.svd, but this will break compatibility.

We should warn/deprecate for a release and then remove both functions for a release. Then we can consider adding torch.Tensor.svd back and mapping it to torch.linalg.svd.

[antocuni]

(sorry guys, I closed the issue by mistake, reopening it)

[pulkin]

You probably do know that numpy returns u, s, vh = svd(A) such that A = u @ s @ vh. I would prefer to have this convention in pytorch because many people use numpy for prototyping.

[antocuni]

while working on PR #45562 we realized that we can't fix torch.svd without breaking existing code. OTOH, the upcoming torch.linalg.svd is completely compatible with numpy.
As discussed on the PR, the current plan is to document the current behavior of torch.svd and tell people to use torch.linalg.svd for new code: this is the commit which tentatively tries to improve the docs, but feel free to suggest a better wording:
abf9baf