From 46a6c0c238dea51fa40796aef6af47b499546c01 Mon Sep 17 00:00:00 2001
From: Joakim Skogholt <onlinereg@skogholt.xyz>
Date: Sun, 14 May 2023 08:49:00 +0200
Subject: [PATCH] Added bidiag struct + decomp function

---
 src/TR.jl | 36 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/src/TR.jl b/src/TR.jl
index 1487c91..147499d 100644
--- a/src/TR.jl
+++ b/src/TR.jl
@@ -22,6 +22,20 @@ struct TRSVD
 end
 
 
+struct TRBidiag
+    W::Matrix{Float64}
+    B::Bidiagonal{Float64, Vector{Float64}}
+    T::Matrix{Float64}
+    y::Vector{Float64}
+    mX::Matrix{Float64}
+    my::Float64
+    regType::String
+    regParam1::Float64
+    regMat::Matrix{Float64}
+    n::Int64
+    p::Int64
+end
+
 
 """
 ### TO DO: ADD FRACTIONAL DERIVATIVE REGULARIZATION ###
@@ -71,6 +85,28 @@ return regMat
 end
 
 
+"""
+    function TRBidiagDecomp(X, y, A=(minimum(size(X))-1), regType="L2", regParam1=0, regParam2=1e-14)
+
+Calculates regularization matrix (using function "RegularizationMatrix"),
+and centres and transforms data matrix according to "X / regMat".
+Output is an object of type "TRBidiag" and is used as input to other TR functions.
+"""
+function TRBidiagDecomp(X, y, A=(minimum(size(X))-1), regType="L2", regParam1=0, regParam2=1e-14)
+
+n, p       = size(X);
+mX         = mean(X, dims=1);    
+X          = X .- mX;
+my         = mean(y);
+y          = vec(y .- my);
+regMat     = regularizationMatrix(X; regType, regParam1, regParam2);
+X          = X / regMat;
+_, W, B, T = bidiag2(X, y, A);
+TRObj      = TRBidiag(W, B, T, y, mX, my, regType, regParam1, regMat, n, p);
+    
+return TRObj
+end
+
 
 """
     function TRSVDDecomp(X, regType="L2", regParam1=0, regParam2=1e-14)