From c4a93d25fa10556c396e6fd05b8ff04be11e487f Mon Sep 17 00:00:00 2001
From: daniel <daniel@kapla.at>
Date: Mon, 5 May 2025 14:34:26 +0200
Subject: [PATCH] add: extended dist.subspace to support list of kronecker
 product components as well as needing less memory, add: if for cond.threshold
 if non-finite, allows to disable regularization, fix: TSIR's mode-wise sample
 covariance scaling factor

---
 tensorPredictors/R/dist_subspace.R      | 53 ++++++++++++++-----------
 tensorPredictors/R/gmlm_tensor_normal.R | 13 +++---
 2 files changed, 37 insertions(+), 29 deletions(-)

diff --git a/tensorPredictors/R/dist_subspace.R b/tensorPredictors/R/dist_subspace.R
index 8fcc9df..935f20b 100644
--- a/tensorPredictors/R/dist_subspace.R
+++ b/tensorPredictors/R/dist_subspace.R
@@ -1,7 +1,7 @@
 #' Subspace distance
 #'
-#' @param A,B Basis matrices as representations of elements of the Grassmann
-#'  manifold.
+#' @param As,Bs Basis matrices or list of Kronecker components of basis matrices
+#' as representations of elements of the Grassmann manifold.
 #' @param is.ortho Boolean to specify if \eqn{A} and \eqn{B} are semi-orthogonal.
 #'   If false, the projection matrices are computed as
 #' \deqn{P_A = A (A' A)^{-1} A'}
@@ -14,36 +14,41 @@
 #'  subspaces of different dimensions" <arXiv:1407.0900>
 #'
 #' @export
-dist.subspace <- function (A, B, is.ortho = FALSE, normalize = FALSE,
+dist.subspace <- function (As, Bs, is.ortho = FALSE, normalize = FALSE,
     tol = sqrt(.Machine$double.eps)
 ) {
-    if (!is.matrix(A)) A <- as.matrix(A)
-    if (!is.matrix(B)) B <- as.matrix(B)
+    As <- if (is.list(As)) Map(as.matrix, As) else list(as.matrix(As))
+    Bs <- if (is.list(Bs)) Map(as.matrix, Bs) else list(as.matrix(Bs))
 
     if (!is.ortho) {
-        qrA <- qr(A, tol)
-        if (qrA$rank < ncol(A)) {
-            A <- qr.Q(qrA)[, abs(diag(qr.R(qrA))) > tol, drop = FALSE]
-        } else {
-            A <- qr.Q(qrA)
-        }
-        qrB <- qr(B, tol)
-        if (qrB$rank < ncol(B)) {
-            B <- qr.Q(qrB)[, abs(diag(qr.R(qrB))) > tol, drop = FALSE]
-        } else {
-            B <- qr.Q(qrB)
-        }
+        As <- Map(function(A) {
+            qrA <- qr(A, tol)
+            if (qrA$rank < ncol(A)) {
+                qr.Q(qrA)[, abs(diag(qr.R(qrA))) > tol, drop = FALSE]
+            } else {
+                qr.Q(qrA)
+            }
+        }, As)
+        Bs <- Map(function(B) {
+            qrB <- qr(B, tol)
+            if (qrB$rank < ncol(B)) {
+                qr.Q(qrB)[, abs(diag(qr.R(qrB))) > tol, drop = FALSE]
+            } else {
+                qr.Q(qrB)
+            }
+        }, Bs)
     }
 
-    PA <- tcrossprod(A, A)
-    PB <- tcrossprod(B, B)
+    rankA <- prod(sapply(As, ncol))
+    rankB <- prod(sapply(Bs, ncol))
 
-    if (normalize) {
-        rankSum <- ncol(A) + ncol(B)
-        c <- 1 / sqrt(max(1, min(rankSum, 2 * nrow(A) - rankSum)))
+    c <- if (normalize) {
+        rankSum <- rankA + rankB
+        sqrt(1 / max(1, min(rankSum, 2 * prod(sapply(As, nrow)) - rankSum)))
     } else {
-        c <- 1
+        1
     }
 
-    c * norm(PA - PB, type = "F")
+    s <- prod(mapply(function(A, B) sum(crossprod(A, B)^2), As, Bs))
+    c * sqrt(max(0, rankA + rankB - 2 * s))
 }
diff --git a/tensorPredictors/R/gmlm_tensor_normal.R b/tensorPredictors/R/gmlm_tensor_normal.R
index 7e7032d..da7c4ac 100644
--- a/tensorPredictors/R/gmlm_tensor_normal.R
+++ b/tensorPredictors/R/gmlm_tensor_normal.R
@@ -105,11 +105,14 @@ gmlm_tensor_normal <- function(X, F, sample.axis = length(dim(X)),
         # Computing `Omega_j`s, the j'th mode presition matrices, in conjunction
         # with regularization of the j'th mode covariance estimate `Sigma_j`
         for (j in seq_along(Sigmas)) {
-            # Compute min and max eigen values
-            min_max <- range(eigen(Sigmas[[j]], TRUE, TRUE)$values)
-            # The condition is approximately `kappa(Sigmas[[j]]) > cond.threshold`
-            if (min_max[2] > cond.threshold * min_max[1]) {
-                Sigmas[[j]] <- Sigmas[[j]] + diag(0.2 * min_max[2], nrow(Sigmas[[j]]))
+            # Regularize Covariances
+            if (is.finite(cond.threshold)) {
+                # Compute min and max eigen values
+                min_max <- range(eigen(Sigmas[[j]], TRUE, TRUE)$values)
+                # The condition is approximately `kappa(Sigmas[[j]]) > cond.threshold`
+                if (min_max[2] > cond.threshold * min_max[1]) {
+                    Sigmas[[j]] <- Sigmas[[j]] + diag(0.2 * min_max[2], nrow(Sigmas[[j]]))
+                }
             }
             # Compute (unconstraint but regularized) Omega_j as covariance inverse
             Omegas[[j]] <- solve(Sigmas[[j]])