Add `abs` matrix function

[araujoms]

I think this is the last matrix function still missing from LinearAlgebra, in the sense that it is both common enough to belong in the standard library and nontrivial enough to deserve a dedicated implementation.

In principle one can define it equally well as abs(A) = sqrt(A'A) or abs(A) = sqrt(A*A'), but the convention is that the former is used; this fits with the convention of writing the polar decomposition as A ≈ U * abs(A) instead of A ≈ abs(A) * U.

The definition is however not a good way to compute it. In the Hermitian case one can simply compute the eigendecomposition and take the absolute value of the eigenvalues; in the general (even non-square) case, probably the best way to compute is via the singular value decomposition: let u, s, v = svd(A), then abs(A) = v' * Diagonal(s) * v. See discussion at https://discourse.julialang.org/t/absolute-value-of-a-matrix/98301

[stevengj]

Another possibility is a Newton-like algorithm for the polar decomposition, e.g. the QDWH algorithm? https://nhigham.com/2015/11/24/faster-svd-via-polar-decomposition/

(At the end, one has $|A| = U^* A$.)

More generally should we just have a function for the polar decomposition?

On the other hand, the SVD-based algorithm is a lot simpler for a first pass.

[stevengj]

Since the polar factorization $A = U |A|$ is a generalization of the polar form $z = e^{i\phi} r$ of a complex number, where abs(z) == r and sign(z) == cis(ϕ) == z / abs(z), should we also have sign(A) == U == A / abs(A)?

[araujoms]

The Newton-like algorithm might very well be faster than the SVD method, one should benchmark it. I'm just worried that it doesn't seem to work for singular matrices.

As for the sign function, I agree that this is the natural definition, but I'm skeptical that users would be able to find it. scipy.linalg.signm is defined differently, and MATLAB doesn't have it at all. The polar decomposition would be nice to have, though.

If we have a sign function then we can compute abs(A) = sign(A) \ A; the other way around, sign(A) = A / abs(A), doesn't work for singular/non-square A.

[stevengj]

doesn't work for singular/non-square A.

A / abs(A) works for non-square A (since the abs(A) definition you are using is always square regardless of A), but I agree that it is problematic for when $A$ does not have full column rank.

[araujoms]

Yeah, the problem is that abs(A) will always be singular when m < n, for a m x n matrix A

One could of course define abs(A) = sqrt(A'A) for m >= n and abs(A) = sqrt(A*A') for m < n, but that's the path to madness.