1

I am making a library with symbolic computations which supports matrices. A matrix may have symbolic entries (eg be over the ring of polynomials of a variable $x$ ie $\mathbb{Q}[x]$).

I have implemented almost all useful linear algebra algorithms and decompositions using exact/symbolic computations (and in some cases fraction-free algorithms) (eg REF, RREF, SNF, INVERSE, PSEUDO-INVERSE, LU, QR, RANK factorisations/decompositions, ROWSPACE, COLUMNSPACE, NULLSPACE etc). The only needed algorithm which I cannot do with exact/symbolic computations is SVD/EVD. I only find numerical algorithms which solve the problem numericaly / approximately and employing "irrational" computations (eg square roots) which are not exact (note: I dont mind if an exact/symbolic algorithm employs square roots, since I can handle these symbolicaly if needed without actually computing square roots).

Any exact/symbolic algorithm for SVD/EVD or any way to compute SVD using one of the decompositions I already have and which are exact?

Note: the library supports arbitrary precision arithmetic (although slow).

An online demo of the library (which is implemented in pure JavaScript and works in browser and node.js) is over here

Nikos M.
  • 2,326

2 Answers2

2

To answer my own question, it seems the answer is:

no exact/symbolic algorithm exists (or is likely to exist) for SVD/EVD

Essentialy the problem is equivalent to the eigenvalue problem:

$$Ax = \lambda x$$

This problem is equivalent to:

$$det(A-\lambda I)=0$$

which is a polynomial equation of $n$th degree in $\lambda$ and according to Abel-Ruffini theorem there is no general exact/closed form solution involving elementary functions and operations (which is what I need in my case, but there is no general closed form solution even if more complex functions are allowed).

Since every monic polynomial can be the characteristic polynomial of a matrix, this means that there is no exact/closed form solution to the eigenvalue problem and no exact/closed form solution to the SVD problem (which can be rephrased as eigenvalue problem) except numerical approximations which unfortunately cannot be used with symbolic entries.

Nikos M.
  • 2,326
  • You write, @Nikos M., that "there is no exact/closed form solution" but, as the first reference claims, "there is no solution in radicals". An exact/closed solution exists, e.g. via algebraic numbers that consists of a polynomial and interval that contains a single root. From this "exact" solution you can extract a numerical approximation using binary search, Newton's method etc. – Dan Feldman Oct 25 '24 at 14:26
  • I was looking for exact closed form symbolic solutions for a software I was writing. But feel free to post another answer. I may be missing something. – Nikos M. Oct 25 '24 at 15:04
0

In the general case, where there is no solution in radicals (as explained here and in many other sites), you can still get exact (closed?) symbolic result as an algebraic number. You can then extract its approximation up to arbitrary desired approximation, or continue with exact (symbolic) computations without adding numerical errors.

I found open code in SymPy, and closed code in Maple. Not sure about Mathematica.

Symbolic eigenvalues and singular values in SymPy:

https://docs.sympy.org/latest/modules/matrices/matrices.html#sympy.matrices.matrixbase.MatrixBase.singular_values

The "algebraic number" is represented via ComplexRootOf.

Even if it cannot be represented as radicals, there is an impressive support for easily extracting approximations up to arbitrary precision from such an algebraic number, via a decimal (float) or a rational numbers. See: https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.rootoftools.ComplexRootOf

In (closed source) Maple, "ComplexRootOf" is similar to "Rootof": https://www.maplesoft.com/support/help/Maple/view.aspx?path=RootOf#:~:text=The%20function%20RootOf%20is%20a,finite%20fields%20(see%20mod).

which is a possible output of symbolic eigenvalue computation: https://www.maplesoft.com/support/help/maple/view.aspx?path=LinearAlgebra%2FEigenvalues

In Mathematica, it seems possible using symbolic input but I could not find an example: https://reference.wolfram.com/language/ref/Eigenvalues.html https://reference.wolfram.com/language/tutorial/ManipulatingEquationsAndInequalities.html#26429

In WolframAlpha I got symbolic solution that I am not sure how to convert, which is why I optimistic regarding (its very related) Mathematica (press the "Exact form" button): https://www.wolframalpha.com/input?i=eigenvalues+%7B%7B4%2C1%2C33%2C4%2C5%7D%2C++%7B4%2C11%2C3%2C4%2C9%7D%2C%7B42%2C1%2C33%2C4%2C8%7D%2C%7B40%2C1%2C3%2C4%2C7%7D%2C%7B4%2C1%2C3%2C4%2C63%7D%7D

Surprisingly, it seems that there is no such support in Sagemath: https://ask.sagemath.org/question/8866/can-i-get-a-exact-solution-for-svd/ https://github.com/sagemath/sage/issues/32171

Dan Feldman
  • 454
  • 2
  • 9
  • As far as I know by studying the source code of sympy, it uses the known exact closed form solutions of polynomials up to 4th degree (which ARE solvable with radicals) to compute exact eigenvalues and singular values, so if the matrix has dimensions >= 5 no such symbolic algorithm exists. But I will check again if now this has changed – Nikos M. Oct 26 '24 at 18:14
  • Indeed singular_values uses eigenvals routine here which explicittly states that it works only for matrices with characteristic polynomials up to 4th degree (which ARE exactly solvable in closed form) – Nikos M. Oct 26 '24 at 18:24
  • Using complexRootOf is another way of saying the solutions of p(x) but without actually having an explicit closed form formula for the solutions which is not what the question asks for. If closed-form formula condition is removed from the question then this method can apply in another sense. – Nikos M. Oct 27 '24 at 09:50
  • You're welcome., 2) there is no common definition for "explicit closed form formula" but algebraic number is the closest candidate in your case. See: https://en.wikipedia.org/wiki/Closed-form_expression. 3) SymPy indeed does not compute eigenvalues for large matrices, but it is easy to use its polynomial solver and compute it yourself as explained in https://stackoverflow.com/questions/65170764/is-it-possible-somehow-to-find-complex-eigenvalues-using-sympy
    • – Dan Feldman Nov 24 '24 at 01:32
  • Besides the closed code that I mentioned, I found out that Sagemath computes symbolically eigenvalues for any matrix (using algebraic numbers, of course).
    • – Dan Feldman Nov 24 '24 at 01:33