SVDComplex -- Compute the SVD decomposition of a chainComplex over RR

Synopsis

Usage:

(h,U)=SVDComplex C or

(h,U)=SVDComplex(C,C')
Inputs:
- C, a chain complex, over RR_{53}
- C', a chain complex, in a lower precision
Optional inputs:
- Strategy => a symbol, default value Projection, Laplacian or Projection for the method used
- Threshold => a real number, default value .0001, the relative threshold used to detect the zero singular values
Outputs:
- h, a hash table, the dimensions of the homology groups HH C
- U, a chain complex map, a map C <- Sigma where the source is the chainComplex of the singular value matrices and U is given by orthogonal matrices

Description

We compute the singular value decomposition either by the iterated Projections or by the Laplacian method. In case the input consists of two chainComplexes we use the iterated Projection method, and identify the stable singular values.

i1 : needsPackage "RandomComplexes"

o1 = RandomComplexes

o1 : Package

i2 : h={1,3,5,2,1}

o2 = {1, 3, 5, 2, 1}

o2 : List

i3 : r={5,11,3,2}

o3 = {5, 11, 3, 2}

o3 : List

i4 : elapsedTime C=randomChainComplex(h,r,Height=>4)
 -- 0.0347177 seconds elapsed

       6       19       19       7       3
o4 = ZZ  <-- ZZ   <-- ZZ   <-- ZZ  <-- ZZ
                                        
     0       1        2        3       4

o4 : ChainComplex

i5 : C.dd^2

           6          19
o5 = 0 : ZZ  <----- ZZ   : 2
                0

           19          7
     1 : ZZ   <----- ZZ  : 3
                 0

           19          3
     2 : ZZ   <----- ZZ  : 4
                 0

o5 : ChainComplexMap

i6 : CR=(C**RR_53)[1]

         6         19         19         7         3
o6 = RR    <-- RR     <-- RR     <-- RR    <-- RR
       53        53         53         53        53
                                                
     -1        0          1          2         3

o6 : ChainComplex

i7 : elapsedTime (h,U)=SVDComplex CR;
 -- 0.0120315 seconds elapsed

i8 : h

o8 = HashTable{-1 => 1}
               0 => 3
               1 => 5
               2 => 2
               3 => 1

o8 : HashTable

i9 : Sigma =source U

         6         19         19         7         3
o9 = RR    <-- RR     <-- RR     <-- RR    <-- RR
       53        53         53         53        53
                                                
     -1        0          1          2         3

o9 : ChainComplex

i10 : Sigma.dd_0

o10 = | 20.7789 0       0       0       0       0 0 0 0 0 0 0 0 0 0 0 0 0 0 |
      | 0       18.3883 0       0       0       0 0 0 0 0 0 0 0 0 0 0 0 0 0 |
      | 0       0       9.51991 0       0       0 0 0 0 0 0 0 0 0 0 0 0 0 0 |
      | 0       0       0       7.19109 0       0 0 0 0 0 0 0 0 0 0 0 0 0 0 |
      | 0       0       0       0       2.40088 0 0 0 0 0 0 0 0 0 0 0 0 0 0 |
      | 0       0       0       0       0       0 0 0 0 0 0 0 0 0 0 0 0 0 0 |

                 6          19
o10 : Matrix RR    <--- RR
               53         53

i11 : errors=apply(toList(min CR+1..max CR),ell->CR.dd_ell-U_(ell-1)*Sigma.dd_ell*transpose U_ell);

i12 : maximalEntry chainComplex errors

o12 = {6.21725e-15, 3.2685e-13, 1.38556e-13, 3.55271e-15}

o12 : List

i13 : elapsedTime (hL,U)=SVDComplex(CR,Strategy=>Laplacian);
 -- 0.0270332 seconds elapsed

i14 : hL === h

o14 = true

i15 : SigmaL =source U;

i16 : for i from min CR+1 to max CR list maximalEntry(SigmaL.dd_i -Sigma.dd_i)

o16 = {1.77636e-14, 2.27374e-13, 5.68434e-14, 1.42109e-14}

o16 : List

i17 : errors=apply(toList(min C+1..max C),ell->CR.dd_ell-U_(ell-1)*SigmaL.dd_ell*transpose U_ell);

i18 : maximalEntry chainComplex errors

o18 = {2.18492e-13, 8.83738e-14, 2.66009e-13, -infinity}

o18 : List

The optional argument

Caveat

The algorithm might fail if the condition numbers of the differential are too bad

Ways to use `SVDComplex` :

"SVDComplex(ChainComplex)"
"SVDComplex(ChainComplex,ChainComplex)"

For the programmer