add FastCholesky for any OPs

X'W-WX == GJG', rank(G) = 2. This allows W to be Cholesky-factored in O(n^2).

Author: MikaelSlevinsky

Project.toml

@@ -1,6 +1,6 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.16.3"
+version = "0.16.4"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/SymmetricToeplitzPlusHankel.jl

@@ -38,23 +38,23 @@ bandwidths(A::SymmetricBandedToeplitzPlusHankel{T}) where T = (A.b, A.b)
 
 #
 # X'W-W*X = G*J*G'
-# This returns G and J, where J = [0 I; -I 0], respecting the skew-symmetry of the right-hand side.
+# This returns G, where J = [0 1; -1 0], respecting the skew-symmetry of the right-hand side.
 #
 function compute_skew_generators(A::SymmetricToeplitzPlusHankel{T}) where T
     v = A.v
     n = size(A, 1)
-    J = [zero(T) one(T); -one(T) zero(T)]
     G = zeros(T, n, 2)
     G[n, 1] = one(T)
-    u2 = reverse(v[2:n+1])
-    u2[1:n-1] .+= v[n+1:2n-1]
-    G[:, 2] .= -u2
-    G, J
+    @inbounds @simd for j in 1:n-1
+        G[j, 2] = -(v[n+2-j] + v[n+j])
+    end
+    G[n, 2] = -v[2]
+    G
 end
 
 function cholesky(A::SymmetricToeplitzPlusHankel{T}) where T
     n = size(A, 1)
-    G, J = compute_skew_generators(A)
+    G = compute_skew_generators(A)
     L = zeros(T, n, n)
     c = A[:, 1]
     ĉ = zeros(T, n)
@@ -74,20 +74,29 @@ function STPHcholesky!(L::Matrix{T}, G, c, ĉ, l, v, row1, n) where T
         for j in 1:n-k+1
             v[j] = G[j, 1]*G[1, 2] - G[j, 2]*G[1, 1]
         end
-        X21 = k == 1 ? T(2) : T(1)
-        ĉ[1] = (c[2] - v[1])/X21
-        for j in 2:n-k
-            ĉ[j] = (c[j+1] + c[j-1] + c[1]*row1[j] - row1[1]*c[j] - v[j])/X21
+        if k == 1
+            for j in 2:n-k
+                ĉ[j] = (c[j+1] + c[j-1] + c[1]*row1[j] - row1[1]*c[j] - v[j])/2
+            end
+            ĉ[n-k+1] = (c[n-k] + c[1]*row1[n-k+1] - row1[1]*c[n-k+1] - v[n-k+1])/2
+            cst = 2/d
+            for j in 1:n-k
+                row1[j] = -cst*l[j+1]
+            end
+        else
+            for j in 2:n-k
+                ĉ[j] = c[j+1] + c[j-1] + c[1]*row1[j] - row1[1]*c[j] - v[j]
+            end
+            ĉ[n-k+1] = c[n-k] + c[1]*row1[n-k+1] - row1[1]*c[n-k+1] - v[n-k+1]
+            cst = 1/d
+            for j in 1:n-k
+                row1[j] = -cst*l[j+1]
+            end
         end
-        ĉ[n-k+1] = (c[n-k] + c[1]*row1[n-k+1] - row1[1]*c[n-k+1] - v[n-k+1])/X21
         cst = c[2]/d
         for j in 1:n-k
             c[j] = ĉ[j+1] - cst*l[j+1]
         end
-        cst = X21/d
-        for j in 1:n-k
-            row1[j] = -cst*l[j+1]
-        end
         gd1 = G[1, 1]/d
         gd2 = G[1, 2]/d
         for j in 1:n-k
@@ -101,37 +110,106 @@ end
 function cholesky(A::SymmetricBandedToeplitzPlusHankel{T}) where T
     n = size(A, 1)
     b = A.b
-    R = BandedMatrix{T}(undef, (n, n), (0, bandwidth(A, 2)))
+    L = BandedMatrix{T}(undef, (n, n), (bandwidth(A, 1), 0))
     c = A[1:b+2, 1]
     ĉ = zeros(T, b+3)
     l = zeros(T, b+3)
     row1 = zeros(T, b+2)
-    SBTPHcholesky!(R, c, ĉ, l, row1, n, b)
-    return Cholesky(R, 'U', 0)
+    SBTPHcholesky!(L, c, ĉ, l, row1, n, b)
+    return Cholesky(L, 'L', 0)
 end
 
-function SBTPHcholesky!(R::BandedMatrix{T}, c, ĉ, l, row1, n, b) where T
+function SBTPHcholesky!(L::BandedMatrix{T}, c, ĉ, l, row1, n, b) where T
     @inbounds @simd for k in 1:n
         d = sqrt(c[1])
         for j in 1:b+1
             l[j] = c[j]/d
         end
         for j in 1:min(n-k+1, b+1)
-            R[k, j+k-1] = l[j]
+            L[j+k-1, k] = l[j]
         end
-        X21 = k == 1 ? T(2) : T(1)
-        ĉ[1] = c[2]/X21
-        for j in 2:b+1
-            ĉ[j] = (c[j+1] + c[j-1] + c[1]*row1[j] - row1[1]*c[j])/X21
+        if k == 1
+            for j in 2:b+1
+                ĉ[j] = (c[j+1] + c[j-1] + c[1]*row1[j] - row1[1]*c[j])/2
+            end
+            ĉ[b+2] = (c[b+1] + c[1]*row1[b+2] - row1[1]*c[b+2])/2
+            cst = 2/d
+            for j in 1:b+2
+                row1[j] = -cst*l[j+1]
+            end
+        else
+            for j in 2:b+1
+                ĉ[j] = (c[j+1] + c[j-1] + c[1]*row1[j] - row1[1]*c[j])
+            end
+            ĉ[b+2] = (c[b+1] + c[1]*row1[b+2] - row1[1]*c[b+2])
+            cst = 1/d
+            for j in 1:b+2
+                row1[j] = -cst*l[j+1]
+            end
         end
-        ĉ[b+2] = (c[b+1] + c[1]*row1[b+2] - row1[1]*c[b+2])/X21
         cst = c[2]/d
         for j in 1:b+2
             c[j] = ĉ[j+1] - cst*l[j+1]
         end
-        cst = X21/d
-        for j in 1:b+2
-            row1[j] = -cst*l[j+1]
+    end
+end
+
+
+
+#
+# X'W-W*X = G*J*G'
+# This returns G, where J = [0 1; -1 0], respecting the skew-symmetry of the right-hand side.
+#
+function compute_skew_generators(W::Symmetric{T, <: AbstractMatrix{T}}, X::Tridiagonal{T, Vector{T}}) where T
+    @assert size(W) == size(X)
+    m, n = size(W)
+    G = zeros(T, n, 2)
+    G[n, 1] = one(T)
+    G[:, 2] .= W[n-1, :]*X[n-1, n] - X'W[:, n]
+    return G
+end
+
+function fastcholesky(W::Symmetric{T, <: AbstractMatrix{T}}, X::Tridiagonal{T, Vector{T}}) where T
+    n = size(W, 1)
+    G = compute_skew_generators(W, X)
+    L = zeros(T, n, n)
+    c = W[:, 1]
+    ĉ = zeros(T, n)
+    l = zeros(T, n)
+    v = zeros(T, n)
+    row1 = zeros(T, n)
+    fastcholesky!(L, X, G, c, ĉ, l, v, row1, n)
+    return Cholesky(L, 'L', 0)
+end
+
+
+function fastcholesky!(L::Matrix{T}, X::Tridiagonal{T, Vector{T}}, G, c, ĉ, l, v, row1, n) where T
+    @inbounds @simd for k in 1:n-1
+        d = sqrt(c[k])
+        for j in k:n
+            L[j, k] = l[j] = c[j]/d
+        end
+        for j in k:n
+            v[j] = G[j, 1]*G[k, 2] - G[j, 2]*G[k, 1]
+        end
+        for j in k+1:n-1
+            ĉ[j] = (X[j-1, j]*c[j-1] + (X[j, j]-X[k, k])*c[j] + X[j+1, j]*c[j+1] + c[k]*row1[j] - row1[k]*c[j] - v[j])/X[k+1, k]
+        end
+        ĉ[n] = (X[n-1, n]*c[n-1] + (X[n, n]-X[k, k])*c[n] + c[k]*row1[n] - row1[k]*c[n] - v[n])/X[k+1, k]
+        cst = X[k+1, k]/d
+        for j in k+1:n
+            row1[j] = -cst*l[j]
+        end
+        cst = c[k+1]/d
+        for j in k:n
+            c[j] = ĉ[j] - cst*l[j]
+        end
+        gd1 = G[k, 1]/d
+        gd2 = G[k, 2]/d
+        for j in k:n
+            G[j, 1] -= l[j]*gd1
+            G[j, 2] -= l[j]*gd2
         end
     end
+    L[n, n] = sqrt(c[n])
 end

test/symmetrictoeplitzplushankeltests.jl

@@ -37,3 +37,15 @@ end
         @test R ≈ U
     end
 end
+
+@testset "Fast Cholesky" begin
+    n = 128
+    for T in (Float32, Float64, BigFloat)
+        R = plan_leg2cheb(T, n; normcheb=true)*I
+        X = Tridiagonal([T(n)/(2n-1) for n in 1:n-1], zeros(T, n), [T(n)/(2n+1) for n in 1:n-1]) # Legendre X
+        W = Symmetric(R'R)
+        F = FastTransforms.fastcholesky(W, X)
+        @test F.L*F.L' ≈ W
+        @test F.U ≈ R
+    end
+end