From a071fbdd360211d66c19f08528fe51794e4c6b80 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Fri, 1 Mar 2024 08:18:45 -0500
Subject: [PATCH] Recursive version that tries to choose a good order

---
 src/gcp-opt.jl | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index e3a5cb2..7fc5412 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -248,12 +248,14 @@ function khatrirao(A::Vararg{T,N}) where {T<:AbstractMatrix,N}
         return A[1]
     end
 
-    # General case: N > 1
-    r = (only ∘ unique)(size.(A, 2))
-    K = similar(A[1], prod(size.(A, 1)), r)
-    for j in 1:r
-        temp = (N == 2) ? view(A[1], :, j) : kron([view(A[i], :, j) for i in 1:N-1]...)
-        kron!(view(K, :, j), temp, view(A[N], :, j))
+    # Base case: N = 2
+    if N == 2
+        r = (only ∘ unique)(size.(A, 2))
+        return reshape(reshape(A[1], :, 1, r) .* reshape(A[2], 1, :, r), :, r)
     end
-    return K
+
+    # Recursive case: N > 2
+    I, r = size.(A, 1), (only ∘ unique)(size.(A, 2))
+    n = argmin(n -> I[n] * I[n+1], 1:N-1)
+    return khatrirao(A[1:n-1]..., khatrirao(A[n], A[n+1]), A[n+2:end]...)
 end