Singular value decomposition

[sylee957]

I have a draft to compute the singular value decomposition with semi-orthogonal U and V, which shows rank-revealing property and which can reduce down a lot of redundant computations I would want to implement this if there are further ideas to reduce the expression blowup issues. (Like skipping normalization)

from sympy import *

def singular_value_decomposition(A):
    AH = A.H
    m, n = A.shape
    if m >= n:
        V, S = (AH * A).diagonalize()

        S_diag = [S[i, i] for i in range(S.rows)]
        ranked = []
        for i, x in enumerate(S.diagonal()):
            if not x.is_zero:
                ranked.append(i)

        V = V[:, ranked]
        S = Matrix.diag([sqrt(S[i, i]) for i in range(S.rows) if i in ranked])

        V, _ = V.QRdecomposition()
        U = A * V * S.inv()
    else:
        U, S = (A * AH).diagonalize()

        S_diag = [S[i, i] for i in range(S.rows)]
        ranked = []
        for i, x in enumerate(S.diagonal()):
            if not x.is_zero:
                ranked.append(i)

        U = U[:, ranked]
        S = Matrix.diag([sqrt(S[i, i]) for i in range(S.rows) if i in ranked])

        U, _ = U.QRdecomposition()
        V = AH * U * S.inv()

    return U, S, V

[safiya03]

Hi @sylee957 I would want to work on this. Can you please help me on how to start.

[ks147]

@sylee957 @oscarbenjamin I was thinking of implementing SVD using polar decomposition. I only saw a function to get the singular values but not for the decomposition. I also didn't see Polar decomposition.

[sylee957]

How does the polar decomposition look like and how should your method be different than this? After I briefly reviewed some materials (e.g. https://web.stanford.edu/~gavish/documents/SVD_ans_you.pdf), if there is no novel way of computing polar decomposition, I think that the singular value decomposition itself can be situated as a step toward polar decomposition because the polar decomposition have the step of evaluating the singular value matrix.

[ks147]

@sylee957, I implemented SVD using Polar Decomposition in one of my courses. I just went through the course material again and realized that it only works for nonsingular matrices. We can find the polar decomposition of a nonsingular matrix A as follows: A = RW, then R * 2 = AA' R = (AA') * (1/2) and W = R.inv() A using this polar decomposition, we can find SVD, but this only works for nonsingular matrices.

You are right; the only way to find SVD for any matrix is as suggested by you originally.

[ks147]

We can easily get Polar decomposition if we already have SVD of a matrix

[ks147]

@sylee957 @oscarbenjamin I think the reason for expression blow-up is the function eigenvals(). When we call the diagonalize function in your code, it makes use of eigenvalues to get diagonalization. and if you take a matrix that doesn't have rational eigenvalues you'll see a blow-up in the Diagonal matrix. The eigenvalues are roots of the characteristic polynomial so they will be represented as a complicated expression if we keep everything in symbolic form obviously.

[sylee957]

We have improved eigenvector computation in #20614, so it would give at least decidable result

[ks147]

My master branch is up to date so I tried running SVD with the given 3 x 3 matrix and I get the following result

from sympy import Matrix A = Matrix([[12, -51, 4], [6, 167, -68], [-4, 24, -41]]) A.SVD() (Matrix([ [ sqrt(6)(4 - 17(-47045644636sqrt(3) + 960750sqrt(10140407) - 143467I(385620038 + 7875sqrt(10140407)I)(2/3) - 47045644636I - 960750sqrt(30421221)I + 102sqrt(10140407)(385620038 + 7875sqrt(10140407)*I)(1/3) + 3sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 143467sqrt(3)(385620038 + 7875sqrt(10140407)I)*(2/3) + 3sqrt(30421221)I(385620038 + 7875sqrt(10140407)I)(2/3) + 209459878I(385620038 + 7875sqrt(10140407)I)(1/3))/(72(-7875sqrt(10140407) + 385620038sqrt(3) + 7875sqrt(30421221)I + 385620038I)) + (3(-34339496 + 24419(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 + sqrt(3)I)2(385620038 + 7875sqrt(10140407)I) + 16(2122156 + (-1534 + (1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)I)(1/3) + 155462782092(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)(2/3))/(47250(1 + sqrt(3)I)2(385620038 + 7875sqrt(10140407)I)))/(7sqrt(1534 - sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3) - (385620038 + 7875sqrt(10140407)I)(1/3) - 2122156/((385620038 + 7875sqrt(10140407)I)(1/3) + sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3)))), sqrt(6)(4 - 17(6(1 - sqrt(3)I)3(886078 + 517(-1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(385620038 + 7875sqrt(10140407)I) + (-1 + sqrt(3)I)(2122156 + (-1534 + (1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)(1/3))(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)I)(1/3) - 6582927912(1 - sqrt(3)I)2(385620038 + 7875sqrt(10140407)*I)(2/3))/(756000(1 - sqrt(3)I)3(385620038 + 7875sqrt(10140407)I)) + (155462782092(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(2/3) + 16(2122156 + (-1534 + (1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)*(1/3))(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)*I)(1/3) + 3(-34339496 + 24419(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 - sqrt(3)I)*2(385620038 + 7875sqrt(10140407)I))/(47250(1 - sqrt(3)I)2(385620038 + 7875sqrt(10140407)I)))/(7sqrt(1534 + 2122156/(-(385620038 + 7875sqrt(10140407)I)(1/3) + sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3)) - (385620038 + 7875sqrt(10140407)I)(1/3) + sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3))), sqrt(3)(4 - 17(-960750sqrt(10140407) - 3sqrt(10140407)(385620038 + 7875sqrt(10140407)*I)(2/3) + 51sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(1/3) + 47045644636I + 143467I(385620038 + 7875sqrt(10140407)*I)(2/3) + 104729939I(385620038 + 7875sqrt(10140407)I)(1/3))/(72(7875sqrt(10140407) - 385620038I)) + (-3811500sqrt(10140407) - 183907I(385620038 + 7875sqrt(10140407)I)(2/3) - 119222294I(385620038 + 7875sqrt(10140407)I)(1/3) - 96sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 72507sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(1/3) + 186640098392I)/(9(-7875sqrt(10140407) + 385620038I)))/(7sqrt(767 + 530539/(385620038 + 7875sqrt(10140407)*I)(1/3) + (385620038 + 7875sqrt(10140407)I)(1/3)))], [sqrt(6)(-68 + 167(-47045644636sqrt(3) + 960750sqrt(10140407) - 143467I(385620038 + 7875sqrt(10140407)I)(2/3) - 47045644636I - 960750sqrt(30421221)I + 102sqrt(10140407)(385620038 + 7875sqrt(10140407)*I)(1/3) + 3sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 143467sqrt(3)(385620038 + 7875sqrt(10140407)I)*(2/3) + 3sqrt(30421221)I(385620038 + 7875sqrt(10140407)I)(2/3) + 209459878I(385620038 + 7875sqrt(10140407)I)(1/3))/(216(-7875sqrt(10140407) + 385620038sqrt(3) + 7875sqrt(30421221)I + 385620038I)) + (3(-34339496 + 24419(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 + sqrt(3)I)2(385620038 + 7875sqrt(10140407)I) + 16(2122156 + (-1534 + (1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)I)(1/3) + 155462782092(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)(2/3))/(94500(1 + sqrt(3)I)2(385620038 + 7875sqrt(10140407)I)))/(7sqrt(1534 - sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3) - (385620038 + 7875sqrt(10140407)I)(1/3) - 2122156/((385620038 + 7875sqrt(10140407)I)(1/3) + sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3)))), sqrt(6)(-68 + 167(6(1 - sqrt(3)I)3(886078 + 517(-1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(385620038 + 7875sqrt(10140407)I) + (-1 + sqrt(3)I)(2122156 + (-1534 + (1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)(1/3))(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)I)(1/3) - 6582927912(1 - sqrt(3)I)2(385620038 + 7875sqrt(10140407)*I)(2/3))/(2268000(1 - sqrt(3)I)3(385620038 + 7875sqrt(10140407)I)) + (155462782092(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(2/3) + 16(2122156 + (-1534 + (1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)*(1/3))(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)*I)(1/3) + 3(-34339496 + 24419(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 - sqrt(3)I)*2(385620038 + 7875sqrt(10140407)I))/(94500(1 - sqrt(3)I)2(385620038 + 7875sqrt(10140407)I)))/(7sqrt(1534 + 2122156/(-(385620038 + 7875sqrt(10140407)I)(1/3) + sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3)) - (385620038 + 7875sqrt(10140407)I)(1/3) + sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3))), sqrt(3)(-68 + (-3811500sqrt(10140407) - 183907I(385620038 + 7875sqrt(10140407)I)(2/3) - 119222294I(385620038 + 7875sqrt(10140407)I)(1/3) - 96sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 72507sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(1/3) + 186640098392I)/(18(-7875sqrt(10140407) + 385620038I)) + 167(-960750sqrt(10140407) - 3sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 51sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(1/3) + 47045644636I + 143467I(385620038 + 7875sqrt(10140407)*I)(2/3) + 104729939I(385620038 + 7875sqrt(10140407)I)(1/3))/(216(7875sqrt(10140407) - 385620038I)))/(7sqrt(767 + 530539/(385620038 + 7875sqrt(10140407)I)(1/3) + (385620038 + 7875sqrt(10140407)I)(1/3)))], [ sqrt(6)(-41 + (-47045644636sqrt(3) + 960750sqrt(10140407) - 143467I(385620038 + 7875sqrt(10140407)*I)(2/3) - 47045644636I - 960750sqrt(30421221)I + 102sqrt(10140407)(385620038 + 7875sqrt(10140407)*I)(1/3) + 3sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 143467sqrt(3)(385620038 + 7875sqrt(10140407)I)*(2/3) + 3sqrt(30421221)I(385620038 + 7875sqrt(10140407)I)(2/3) + 209459878I(385620038 + 7875sqrt(10140407)I)(1/3))/(9(-7875sqrt(10140407) + 385620038sqrt(3) + 7875sqrt(30421221)I + 385620038I)) - (3(-34339496 + 24419(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 + sqrt(3)I)2(385620038 + 7875sqrt(10140407)I) + 16(2122156 + (-1534 + (1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)I)(1/3) + 155462782092(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)(2/3))/(141750(1 + sqrt(3)I)2(385620038 + 7875sqrt(10140407)I)))/(7sqrt(1534 - sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3) - (385620038 + 7875sqrt(10140407)I)(1/3) - 2122156/((385620038 + 7875sqrt(10140407)I)(1/3) + sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3)))), sqrt(6)(-41 + (6(1 - sqrt(3)*I)3(886078 + 517(-1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(385620038 + 7875sqrt(10140407)I) + (-1 + sqrt(3)I)(2122156 + (-1534 + (1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)*(1/3))(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)*I)(1/3) - 6582927912(1 - sqrt(3)I)2(385620038 + 7875sqrt(10140407)*I)(2/3))/(94500(1 - sqrt(3)I)3(385620038 + 7875sqrt(10140407)I)) - (155462782092(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(2/3) + 16(2122156 + (-1534 + (1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)(1/3))(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)I)(1/3) + 3(-34339496 + 24419(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 - sqrt(3)I)2(385620038 + 7875sqrt(10140407)I))/(141750(1 - sqrt(3)*I)2(385620038 + 7875sqrt(10140407)I)))/(7sqrt(1534 + 2122156/(-(385620038 + 7875sqrt(10140407)I)(1/3) + sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3)) - (385620038 + 7875sqrt(10140407)I)(1/3) + sqrt(3)I(385620038 + 7875sqrt(10140407)I)(1/3))), sqrt(3)(-41 - (-3811500sqrt(10140407) - 183907I(385620038 + 7875sqrt(10140407)I)(2/3) - 119222294I(385620038 + 7875sqrt(10140407)I)(1/3) - 96sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 72507sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(1/3) + 186640098392I)/(27(-7875sqrt(10140407) + 385620038I)) + (-960750sqrt(10140407) - 3sqrt(10140407)(385620038 + 7875sqrt(10140407)*I)(2/3) + 51sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(1/3) + 47045644636I + 143467I(385620038 + 7875sqrt(10140407)*I)(2/3) + 104729939I(385620038 + 7875sqrt(10140407)I)(1/3))/(9(7875sqrt(10140407) - 385620038I)))/(7sqrt(767 + 530539/(385620038 + 7875sqrt(10140407)I)(1/3) + (385620038 + 7875sqrt(10140407)I)(1/3)))]]), Matrix([ [sqrt(37583/3 + 49(-1/2 - sqrt(3)I/2)(385620038/27 + 875sqrt(10140407)*I/3)*(1/3) + 25996411/(9(-1/2 - sqrt(3)I/2)(385620038/27 + 875sqrt(10140407)I/3)(1/3))), 0, 0], [ 0, sqrt(37583/3 + 25996411/(9(-1/2 + sqrt(3)I/2)(385620038/27 + 875sqrt(10140407)*I/3)(1/3)) + 49(-1/2 + sqrt(3)I/2)(385620038/27 + 875sqrt(10140407)*I/3)*(1/3)), 0], [ 0, 0, sqrt(37583/3 + 25996411/(9(385620038/27 + 875sqrt(10140407)I/3)(1/3)) + 49(385620038/27 + 875sqrt(10140407)*I/3)(1/3))]]), Matrix([ [ (3(-34339496 + 24419(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 + sqrt(3)I)2(385620038 + 7875sqrt(10140407)I) + 16(2122156 + (-1534 + (1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)I)(1/3) + 155462782092(1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)(2/3))/(567000(1 + sqrt(3)I)2(385620038 + 7875sqrt(10140407)I)), (155462782092(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)*(2/3) + 16(2122156 + (-1534 + (1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)I)(1/3) + 3(-34339496 + 24419(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(1 - sqrt(3)I)2(385620038 + 7875sqrt(10140407)I))/(567000(1 - sqrt(3)*I)2(385620038 + 7875sqrt(10140407)I)), (-3811500sqrt(10140407) - 183907I(385620038 + 7875sqrt(10140407)I)(2/3) - 119222294I(385620038 + 7875sqrt(10140407)I)(1/3) - 96sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 72507sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(1/3) + 186640098392I)/(108(-7875sqrt(10140407) + 385620038I))], [(-47045644636sqrt(3) + 960750sqrt(10140407) - 143467I(385620038 + 7875sqrt(10140407)I)(2/3) - 47045644636I - 960750sqrt(30421221)I + 102sqrt(10140407)(385620038 + 7875sqrt(10140407)*I)(1/3) + 3sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 143467sqrt(3)(385620038 + 7875sqrt(10140407)I)(2/3) + 3sqrt(30421221)I(385620038 + 7875sqrt(10140407)*I)(2/3) + 209459878I(385620038 + 7875sqrt(10140407)I)(1/3))/(216(-7875sqrt(10140407) + 385620038sqrt(3) + 7875sqrt(30421221)I + 385620038I)), (6(1 - sqrt(3)I)3(886078 + 517(-1 + sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))(385620038 + 7875sqrt(10140407)I) + (-1 + sqrt(3)I)(2122156 + (-1534 + (1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)*I)(1/3))(1 - sqrt(3)I)(385620038 + 7875sqrt(10140407)I)(1/3))2(385620038 + 7875sqrt(10140407)I)(1/3) - 6582927912(1 - sqrt(3)I)2(385620038 + 7875sqrt(10140407)*I)(2/3))/(2268000(1 - sqrt(3)I)3(385620038 + 7875sqrt(10140407)I)), (-960750sqrt(10140407) - 3sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(2/3) + 51sqrt(10140407)(385620038 + 7875sqrt(10140407)I)(1/3) + 47045644636I + 143467I(385620038 + 7875sqrt(10140407)*I)(2/3) + 104729939I(385620038 + 7875sqrt(10140407)I)*(1/3))/(216(7875sqrt(10140407) - 385620038I))], [ 1, 1, 1]]))

[ks147]

__(385620038 + 7875_sqrt(10140407)I))/(141750(1 - sqrt(3)_I)2(385620038 + 7875_sqrt(10140407)_I)))/(7_sqrt(1534 + 2122156/(-(385620038 + 7875_sqrt(10140407)_I)(1/3) + sqrt(3)I(385620038 + 7875_sqrt(10140407)_I)(1/3)) - (385620038 + 7875_sqrt(10140407)_I)(1/3) + sqrt(3)I(385620038 + 7875_sqrt(10140407)_I)(1/3))), sqrt(3)_(-41 - (-3811500_sqrt(10140407) - 183907_I(385620038 + 7875_sqrt(10140407)_I)(2/3) - 119222294I(385620038 + 7875_sqrt(10140407)_I)(1/3) - 96sqrt(10140407)(385620038 + 7875_sqrt(10140407)_I)(2/3) + 72507sqrt(10140407)(385620038 + 7875_sqrt(10140407)_I)(1/3) + 186640098392I)/(27(-7875_sqrt(10140407) + 385620038_I)) + (-960750_sqrt(10140407) - 3sqrt(10140407)(385620038 + 7875_sqrt(10140407)_I)(2/3) + 51sqrt(10140407)(385620038 + 7875_sqrt(10140407)_I)(1/3) + 47045644636_I + 143467I(385620038 + 7875_sqrt(10140407)_I)(2/3) + 104729939I(385620038 + 7875_sqrt(10140407)_I)(1/3))/(9_(7875_sqrt(10140407) - 385620038_I)))/(7_sqrt(767 + 530539/(385620038 + 7875_sqrt(10140407)_I)(1/3) + (385620038 + 7875_sqrt(10140407)_I)(1/3)))]]), Matrix([ [sqrt(37583/3 + 49_(-1/2 - sqrt(3)I/2)(385620038/27 + 875_sqrt(10140407)I/3)**(1/3) + 25996411/(9(-1/2 - sqrt(3)I/2)(385620038/27 + 875_sqrt(10140407)I/3)**(1/3))), 0, 0], [ 0, sqrt(37583/3 + 25996411/(9(-1/2 + sqrt(3)I/2)(385620038/27 + 875_sqrt(10140407)I/3)**(1/3)) + 49(-1/2 + sqrt(3)I/2)(385620038/27 + 875_sqrt(10140407)I/3)**(1/3)), 0], [ 0, 0, sqrt(37583/3 + 25996411/(9(385620038/27 + 875_sqrt(10140407)I/3)**(1/3)) + 49(385620038/27 + 875_sqrt(10140407)I/3)**(1/3))]]), Matrix([ [ (3(-34339496 + 24419_(1 + sqrt(3)I)(385620038 + 7875sqrt(1014

@sylee957 If this is fine with you'll I'll send a PR for SVD and after SVD I can start working on other related things like Polar decomposition