From f869de6a3765fdbb25a6c914463e13fdaeac170d Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Tue, 28 Nov 2023 09:23:08 -0500
Subject: [PATCH 01/40] Fixes #33

---
 src/gcp-opt.jl | 43 +++++++++++++++++++++++++++++--------------
 1 file changed, 29 insertions(+), 14 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 4716326..2a64f44 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -18,6 +18,10 @@ If the LossFunctions.jl package is also loaded,
 Check `GCPDecompositions.LossFunctionsExt.SupportedLosses`
 to see what losses are supported.
 
+Algorithm for computing one mode of MTTKRP is from "Fast Alternating LS Algorithms
+for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al., specifically
+section III-B.
+
 See also: `CPD`, `AbstractLoss`.
 """
 gcp(
@@ -65,7 +69,6 @@ function _gcp(
 ) where {TX,N}
     # T = promote_type(nonmissingtype(TX), Float64)
     T = Float64    # LBFGSB.jl seems to only support Float64
-
     # Compute lower bound from constraints
     lower = maximum(constraint.value for constraint in constraints; init = T(-Inf))
 
@@ -93,7 +96,7 @@ function _gcp(
         M0.U[k] .*= (Xnorm / M0norm)^(1 / N)
     end
     u0 = vcat(vec.(M0.U)...)
-
+    
     # Setup vectorized objective function and gradient
     vec_cutoffs = (0, cumsum(r .* size(X))...)
     vec_ranges = ntuple(k -> vec_cutoffs[k]+1:vec_cutoffs[k+1], Val(N))
@@ -107,7 +110,6 @@ function _gcp(
         gcp_grad_U!(GU, CPD(ones(T, r), U), X, loss)
         return gu
     end
-
     # Run LBFGSB
     lbfgsopts = (; (pn => getproperty(algorithm, pn) for pn in propertynames(algorithm))...)
     u = lbfgsb(f, g!, u0; lb = fill(lower, length(u0)), lbfgsopts...)[2]
@@ -174,21 +176,34 @@ function _gcp(
     return CPD(λ, Tuple(U))
 end
 
-# inefficient but simple
+# Faster MTTKRP
 function mttkrp(X, U, n)
+
     # Dimensions
     N, I, r = length(U), Tuple(size.(U, 1)), (only∘unique)(size.(U, 2))
     (N == ndims(X) && I == size(X)) || throw(DimensionMismatch("`X` and `U` do not have matching dimensions"))
 
-    # Matricized tensor (in mode n)
-    Xn = reshape(permutedims(X, [n; setdiff(1:N, n)]), size(X, n), :)
-
-    # Khatri-Rao product (in mode n)
-    Zn = similar(U[1], prod(I[setdiff(1:N, n)]), r)
+    # See section III-B from "Fast Alternating LS Algorithms for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al.
+    Rn = similar(U[n])
+    Jn = prod(size(X)[1:n])
+    Kn = prod(size(X)[n+1:end])
+    # Compute tensor-vector products right to left (equations 15, 17) for each rank
     for j in 1:r
-        Zn[:, j] = reduce(kron, [view(U[i], :, j) for i in reverse(setdiff(1:N, n))])
+        # Inner tensor-vector products
+        Rn_j = reshape(reshape(X, Jn, Kn) * reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)], init=1), size(X)[1:n]) 
+        
+        # Outer tensor-vector products
+        # Permute so dims to be multiplied are last
+        Rn_j = permutedims(Rn_j, [n; 1:n-1])
+        # Multiply from right to left
+        sz = size(Rn_j)
+        m = length(sz)
+        for k in n-1:-1:1
+            Rn_j = reshape(Rn_j, prod(sz[1:m-1]), sz[m])
+            Rn_j = Rn_j * U[k][:, j]
+            m -= 1
+        end
+        Rn[:, j] = Rn_j  
     end
-
-    # MTTKRP (in mode n)
-    return Xn * Zn
-end
+    return Rn
+end
\ No newline at end of file

From a933580271b24b3c7b70532a7887062dca3e563f Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Tue, 28 Nov 2023 09:53:12 -0500
Subject: [PATCH 02/40] Modify computation outer tensor-vector products

---
 src/gcp-opt.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 2a64f44..b0e9c45 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -199,8 +199,7 @@ function mttkrp(X, U, n)
         sz = size(Rn_j)
         m = length(sz)
         for k in n-1:-1:1
-            Rn_j = reshape(Rn_j, prod(sz[1:m-1]), sz[m])
-            Rn_j = Rn_j * U[k][:, j]
+            Rn_j = reshape(Rn_j, prod(sz[1:m-1]), sz[m]) * U[k][:, j]
             m -= 1
         end
         Rn[:, j] = Rn_j  

From 060173a55830b7f46c0664148bcea90e2be81f60 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Tue, 28 Nov 2023 14:43:38 -0500
Subject: [PATCH 03/40] Do outer tensor-vector products more efficiently

---
 src/gcp-opt.jl | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index b0e9c45..8fcb110 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -190,19 +190,22 @@ function mttkrp(X, U, n)
     # Compute tensor-vector products right to left (equations 15, 17) for each rank
     for j in 1:r
         # Inner tensor-vector products
-        Rn_j = reshape(reshape(X, Jn, Kn) * reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)], init=1), size(X)[1:n]) 
-        
+        inner = reshape(reshape(X, Jn, Kn) * reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)], init=1), size(X)[1:n]) 
         # Outer tensor-vector products
+        Jn_inner = prod(size(inner)[1:n-1])
+        Kn_inner = prod(size(inner)[n:end])
+        Rn[:, j] = transpose(reshape(inner, Jn_inner, Kn_inner)) * reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)], 
+                        init=1)
         # Permute so dims to be multiplied are last
-        Rn_j = permutedims(Rn_j, [n; 1:n-1])
+        #Rn_j = permutedims(Rn_j, [n; 1:n-1])
         # Multiply from right to left
-        sz = size(Rn_j)
-        m = length(sz)
-        for k in n-1:-1:1
-            Rn_j = reshape(Rn_j, prod(sz[1:m-1]), sz[m]) * U[k][:, j]
-            m -= 1
-        end
-        Rn[:, j] = Rn_j  
+        #sz = size(Rn_j)
+        #m = length(sz)
+        #for k in n-1:-1:1
+        #    Rn_j = reshape(Rn_j, prod(sz[1:m-1]), sz[m]) * U[k][:, j]
+        #    m -= 1
+        #end
+        #Rn[:, j] = Rn_j  
     end
     return Rn
 end
\ No newline at end of file

From 08e6323bb6030f3b056011b29f410d4e9f954e16 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Tue, 28 Nov 2023 16:41:15 -0500
Subject: [PATCH 04/40] Clean up

---
 src/gcp-opt.jl | 25 ++++++++-----------------
 1 file changed, 8 insertions(+), 17 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 8fcb110..383ba55 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -18,10 +18,6 @@ If the LossFunctions.jl package is also loaded,
 Check `GCPDecompositions.LossFunctionsExt.SupportedLosses`
 to see what losses are supported.
 
-Algorithm for computing one mode of MTTKRP is from "Fast Alternating LS Algorithms
-for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al., specifically
-section III-B.
-
 See also: `CPD`, `AbstractLoss`.
 """
 gcp(
@@ -176,7 +172,13 @@ function _gcp(
     return CPD(λ, Tuple(U))
 end
 
-# Faster MTTKRP
+"""
+    mttkrp(X, U, n) -> Rn
+    
+    Algorithm for computing one mode of MTTKRP is from "Fast Alternating LS Algorithms
+    for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al., specifically
+    section III-B.
+"""
 function mttkrp(X, U, n)
 
     # Dimensions
@@ -194,18 +196,7 @@ function mttkrp(X, U, n)
         # Outer tensor-vector products
         Jn_inner = prod(size(inner)[1:n-1])
         Kn_inner = prod(size(inner)[n:end])
-        Rn[:, j] = transpose(reshape(inner, Jn_inner, Kn_inner)) * reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)], 
-                        init=1)
-        # Permute so dims to be multiplied are last
-        #Rn_j = permutedims(Rn_j, [n; 1:n-1])
-        # Multiply from right to left
-        #sz = size(Rn_j)
-        #m = length(sz)
-        #for k in n-1:-1:1
-        #    Rn_j = reshape(Rn_j, prod(sz[1:m-1]), sz[m]) * U[k][:, j]
-        #    m -= 1
-        #end
-        #Rn[:, j] = Rn_j  
+        Rn[:, j] = transpose(reshape(inner, Jn_inner, Kn_inner)) * reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)], init=1)
     end
     return Rn
 end
\ No newline at end of file

From 72f88c80e802d49cf67b23f224474e969f1cc8f7 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Mon, 4 Dec 2023 16:33:17 -0500
Subject: [PATCH 05/40] Split into special cases

---
 src/gcp-opt.jl | 28 ++++++++++++++++++++++------
 1 file changed, 22 insertions(+), 6 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 383ba55..cf680cc 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -191,12 +191,28 @@ function mttkrp(X, U, n)
     Kn = prod(size(X)[n+1:end])
     # Compute tensor-vector products right to left (equations 15, 17) for each rank
     for j in 1:r
-        # Inner tensor-vector products
-        inner = reshape(reshape(X, Jn, Kn) * reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)], init=1), size(X)[1:n]) 
-        # Outer tensor-vector products
-        Jn_inner = prod(size(inner)[1:n-1])
-        Kn_inner = prod(size(inner)[n:end])
-        Rn[:, j] = transpose(reshape(inner, Jn_inner, Kn_inner)) * reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)], init=1)
+
+        # Special cases are n = 1 and n = N (n = 1 has no outer tensor-vector products),
+        # n = N has no inner tensor-vector products
+        if n == 1
+            # Just inner tensor-vector products
+            kr_inner = reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)])
+            Rn[:, j] = transpose(reshape(X, Jn, Kn) * kr_inner)
+        elseif n == N
+            # Just outer tensor-vector products
+            kr_outer = reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)])
+            Rn[:, j] = transpose(reshape(X, prod(size(X)[1:n-1]), prod(size(X)[n:end]))) * kr_outer
+        else
+            # Inner tensor-vector products
+            kr_inner = reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)])
+            inner = reshape(reshape(X, Jn, Kn) * kr_inner, size(X)[1:n]) 
+            # Outer tensor-vector products
+            Jn_inner = prod(size(inner)[1:n-1])
+            Kn_inner = prod(size(inner)[n:end])
+            kr_outer = reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)])
+            Rn[:, j] = transpose(reshape(inner, Jn_inner, Kn_inner)) * kr_outer
+        end
+        
     end
     return Rn
 end
\ No newline at end of file

From 17f745ca7c9631dfdf6bba9d9b585f731d77c73f Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Fri, 5 Jan 2024 17:46:38 -0500
Subject: [PATCH 06/40] Nest for loop inside split into special cases

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

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index cf680cc..b72df48 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -190,19 +190,23 @@ function mttkrp(X, U, n)
     Jn = prod(size(X)[1:n])
     Kn = prod(size(X)[n+1:end])
     # Compute tensor-vector products right to left (equations 15, 17) for each rank
-    for j in 1:r
 
-        # Special cases are n = 1 and n = N (n = 1 has no outer tensor-vector products),
-        # n = N has no inner tensor-vector products
-        if n == 1
+    # Special cases are n = 1 and n = N (n = 1 has no outer tensor-vector products),
+    # n = N has no inner tensor-vector products
+    if n == 1
+        for j in 1:r
             # Just inner tensor-vector products
             kr_inner = reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)])
             Rn[:, j] = transpose(reshape(X, Jn, Kn) * kr_inner)
-        elseif n == N
+        end
+    elseif n == N
+        for j in 1:r
             # Just outer tensor-vector products
             kr_outer = reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)])
             Rn[:, j] = transpose(reshape(X, prod(size(X)[1:n-1]), prod(size(X)[n:end]))) * kr_outer
-        else
+        end
+    else
+        for j in 1:r
             # Inner tensor-vector products
             kr_inner = reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)])
             inner = reshape(reshape(X, Jn, Kn) * kr_inner, size(X)[1:n]) 
@@ -212,7 +216,6 @@ function mttkrp(X, U, n)
             kr_outer = reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)])
             Rn[:, j] = transpose(reshape(inner, Jn_inner, Kn_inner)) * kr_outer
         end
-        
     end
     return Rn
 end
\ No newline at end of file

From 8987dbd53e25bf1bfa6e607fb1eb1debc54eda83 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Fri, 5 Jan 2024 17:49:36 -0500
Subject: [PATCH 07/40] Update benchmark manifest to Julia 1.10

---
 benchmark/Manifest.toml | 44 +++++++++++++++++++++++------------------
 1 file changed, 25 insertions(+), 19 deletions(-)

diff --git a/benchmark/Manifest.toml b/benchmark/Manifest.toml
index 4b618f9..1d840f6 100644
--- a/benchmark/Manifest.toml
+++ b/benchmark/Manifest.toml
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.9.3"
+julia_version = "1.10.0"
 manifest_format = "2.0"
 project_hash = "180b24a154cbc6c49a70194d9e5abb3ac49bf9b8"
 
@@ -84,7 +84,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+0"
+version = "1.0.5+1"
 
 [[deps.Contour]]
 git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781"
@@ -265,21 +265,26 @@ version = "0.2.0"
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
-version = "0.6.3"
+version = "0.6.4"
 
 [[deps.LibCURL_jll]]
 deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
-version = "7.84.0+0"
+version = "8.4.0+0"
 
 [[deps.LibGit2]]
-deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
+[[deps.LibGit2_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"]
+uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5"
+version = "1.6.4+0"
+
 [[deps.LibSSH2_jll]]
 deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
-version = "1.10.2+0"
+version = "1.11.0+1"
 
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -326,7 +331,7 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.28.2+0"
+version = "2.28.2+1"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
@@ -339,7 +344,7 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2022.10.11"
+version = "2023.1.10"
 
 [[deps.NaNMath]]
 deps = ["OpenLibm_jll"]
@@ -354,12 +359,12 @@ version = "1.2.0"
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.21+4"
+version = "0.3.23+2"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.1+0"
+version = "0.8.1+2"
 
 [[deps.OpenSpecFun_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
@@ -387,7 +392,7 @@ version = "2.7.2"
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-version = "1.9.2"
+version = "1.10.0"
 
 [[deps.PkgBenchmark]]
 deps = ["BenchmarkTools", "Dates", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Pkg", "Printf", "TerminalLoggers", "UUIDs"]
@@ -432,7 +437,7 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[deps.Random]]
-deps = ["SHA", "Serialization"]
+deps = ["SHA"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[deps.Reexport]]
@@ -477,6 +482,7 @@ version = "1.2.0"
 [[deps.SparseArrays]]
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+version = "1.10.0"
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
@@ -514,7 +520,7 @@ version = "1.4.2"
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.StatsAPI]]
 deps = ["LinearAlgebra"]
@@ -547,9 +553,9 @@ deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
 uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [[deps.SuiteSparse_jll]]
-deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
+deps = ["Artifacts", "Libdl", "libblastrampoline_jll"]
 uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
-version = "5.10.1+6"
+version = "7.2.1+1"
 
 [[deps.TOML]]
 deps = ["Dates"]
@@ -613,19 +619,19 @@ version = "3.6.0"
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.13+0"
+version = "1.2.13+1"
 
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+0"
+version = "5.8.0+1"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
-version = "1.48.0+0"
+version = "1.52.0+1"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
-version = "17.4.0+0"
+version = "17.4.0+2"

From bee7da2d16862c218a02d7b20a5592a291952fca Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Fri, 5 Jan 2024 17:59:32 -0500
Subject: [PATCH 08/40] Amortize the matrix multiply for n=1 case

Precompute (and store) the full khatri-rao product then do one matrix multiply
---
 src/gcp-opt.jl | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index b72df48..f1fe72a 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -194,11 +194,12 @@ function mttkrp(X, U, n)
     # Special cases are n = 1 and n = N (n = 1 has no outer tensor-vector products),
     # n = N has no inner tensor-vector products
     if n == 1
+        # Just inner tensor-vector products
+        kr_inner = similar(U[1], prod(I[2:N]), r)
         for j in 1:r
-            # Just inner tensor-vector products
-            kr_inner = reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)])
-            Rn[:, j] = transpose(reshape(X, Jn, Kn) * kr_inner)
+            kr_inner[:,j] = reduce(kron, [view(U[i], :, j) for i in reverse(2:N)])
         end
+        Rn = reshape(X, size(X, 1), :) * kr_inner
     elseif n == N
         for j in 1:r
             # Just outer tensor-vector products

From 9d1fc8d67d0fc054135e00834d7470ba64410443 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Fri, 5 Jan 2024 19:25:26 -0500
Subject: [PATCH 09/40] Revert "Update benchmark manifest to Julia 1.10"

This reverts commit 8987dbd53e25bf1bfa6e607fb1eb1debc54eda83.
---
 benchmark/Manifest.toml | 44 ++++++++++++++++++-----------------------
 1 file changed, 19 insertions(+), 25 deletions(-)

diff --git a/benchmark/Manifest.toml b/benchmark/Manifest.toml
index 1d840f6..4b618f9 100644
--- a/benchmark/Manifest.toml
+++ b/benchmark/Manifest.toml
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.9.3"
 manifest_format = "2.0"
 project_hash = "180b24a154cbc6c49a70194d9e5abb3ac49bf9b8"
 
@@ -84,7 +84,7 @@ weakdeps = ["Dates", "LinearAlgebra"]
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.0.5+0"
 
 [[deps.Contour]]
 git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781"
@@ -265,26 +265,21 @@ version = "0.2.0"
 [[deps.LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
-version = "0.6.4"
+version = "0.6.3"
 
 [[deps.LibCURL_jll]]
 deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
-version = "8.4.0+0"
+version = "7.84.0+0"
 
 [[deps.LibGit2]]
-deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
-[[deps.LibGit2_jll]]
-deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"]
-uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5"
-version = "1.6.4+0"
-
 [[deps.LibSSH2_jll]]
 deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
-version = "1.11.0+1"
+version = "1.10.2+0"
 
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -331,7 +326,7 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 [[deps.MbedTLS_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.28.2+1"
+version = "2.28.2+0"
 
 [[deps.Missings]]
 deps = ["DataAPI"]
@@ -344,7 +339,7 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2023.1.10"
+version = "2022.10.11"
 
 [[deps.NaNMath]]
 deps = ["OpenLibm_jll"]
@@ -359,12 +354,12 @@ version = "1.2.0"
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.21+4"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
-version = "0.8.1+2"
+version = "0.8.1+0"
 
 [[deps.OpenSpecFun_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
@@ -392,7 +387,7 @@ version = "2.7.2"
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-version = "1.10.0"
+version = "1.9.2"
 
 [[deps.PkgBenchmark]]
 deps = ["BenchmarkTools", "Dates", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Pkg", "Printf", "TerminalLoggers", "UUIDs"]
@@ -437,7 +432,7 @@ deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[deps.Random]]
-deps = ["SHA"]
+deps = ["SHA", "Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [[deps.Reexport]]
@@ -482,7 +477,6 @@ version = "1.2.0"
 [[deps.SparseArrays]]
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-version = "1.10.0"
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
@@ -520,7 +514,7 @@ version = "1.4.2"
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-version = "1.10.0"
+version = "1.9.0"
 
 [[deps.StatsAPI]]
 deps = ["LinearAlgebra"]
@@ -553,9 +547,9 @@ deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
 uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [[deps.SuiteSparse_jll]]
-deps = ["Artifacts", "Libdl", "libblastrampoline_jll"]
+deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
 uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
-version = "7.2.1+1"
+version = "5.10.1+6"
 
 [[deps.TOML]]
 deps = ["Dates"]
@@ -619,19 +613,19 @@ version = "3.6.0"
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.13+1"
+version = "1.2.13+0"
 
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
-version = "5.8.0+1"
+version = "5.8.0+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
-version = "1.52.0+1"
+version = "1.48.0+0"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
-version = "17.4.0+2"
+version = "17.4.0+0"

From 99698b517092a8fcea3ad91e5777ecc7f1d7818f Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Fri, 5 Jan 2024 19:37:09 -0500
Subject: [PATCH 10/40] Factor out khatri-rao product

---
 src/gcp-opt.jl | 22 +++++++++++++++++-----
 1 file changed, 17 insertions(+), 5 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index f1fe72a..2d6419c 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -195,10 +195,7 @@ function mttkrp(X, U, n)
     # n = N has no inner tensor-vector products
     if n == 1
         # Just inner tensor-vector products
-        kr_inner = similar(U[1], prod(I[2:N]), r)
-        for j in 1:r
-            kr_inner[:,j] = reduce(kron, [view(U[i], :, j) for i in reverse(2:N)])
-        end
+        kr_inner = khatrirao(U[reverse(2:N)]...)
         Rn = reshape(X, size(X, 1), :) * kr_inner
     elseif n == N
         for j in 1:r
@@ -219,4 +216,19 @@ function mttkrp(X, U, n)
         end
     end
     return Rn
-end
\ No newline at end of file
+end
+
+"""
+    khatrirao(A1, A2, ...)
+    
+    Computes the Khatri-Rao product (i.e., column-wise Kronecker product)
+    of the matrices `A1`, `A2`, etc.
+"""
+function khatrirao(A::Vararg{T,N}) where {T<:Matrix,N}
+    r = (only∘unique)(size.(A,2))
+    K = similar(A[1], prod(size.(A,1)), r)
+    for j in 1:r
+        K[:, j] = reduce(kron, [view(A[i], :, j) for i in 1:N])
+    end
+    return K
+end

From aec931cdccbf27760bcef7a8abb20070838e3874 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Fri, 5 Jan 2024 19:52:53 -0500
Subject: [PATCH 11/40] Use in-place `kron!` to (slightly) reduce allocations

---
 src/gcp-opt.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 2d6419c..d84f614 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -228,7 +228,12 @@ function khatrirao(A::Vararg{T,N}) where {T<:Matrix,N}
     r = (only∘unique)(size.(A,2))
     K = similar(A[1], prod(size.(A,1)), r)
     for j in 1:r
-        K[:, j] = reduce(kron, [view(A[i], :, j) for i in 1:N])
+        # Quick way to avoid making two copies of `K`:
+        # 1. use `reduce` to compute all but the last `kron`
+        # 2. use `kron!` for the last one
+        # Can/should optimize later: https://github.com/dahong67/GCPDecompositions.jl/issues/34
+        temp = reduce(kron, [view(A[i], :, j) for i in 1:N-1])
+        kron!(view(K, :, j), temp, view(A[N], :, j))
     end
     return K
 end

From ad785debe22edafe24bbfca71081f203b93024d7 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Mon, 8 Jan 2024 10:26:11 -0500
Subject: [PATCH 12/40] Use khatri-rao function, combine matrix-vector
 multiplies for n=N case

---
 src/gcp-opt.jl | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index d84f614..e4b9283 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -198,11 +198,9 @@ function mttkrp(X, U, n)
         kr_inner = khatrirao(U[reverse(2:N)]...)
         Rn = reshape(X, size(X, 1), :) * kr_inner
     elseif n == N
-        for j in 1:r
-            # Just outer tensor-vector products
-            kr_outer = reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)])
-            Rn[:, j] = transpose(reshape(X, prod(size(X)[1:n-1]), prod(size(X)[n:end]))) * kr_outer
-        end
+        # Just outer tensor-vector products
+        kr_outer = khatrirao(U[reverse(1:N-1)]...)
+        Rn = transpose(reshape(X, prod(size(X)[1:N-1]), size(X)[N])) * kr_outer
     else
         for j in 1:r
             # Inner tensor-vector products

From d5ffc8585274fe89d0077e3454b6fe758090ccd4 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Mon, 8 Jan 2024 12:08:02 -0500
Subject: [PATCH 13/40] Use in-place matrix-matrix multiplies

---
 src/gcp-opt.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index e4b9283..37b2daa 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -196,11 +196,11 @@ function mttkrp(X, U, n)
     if n == 1
         # Just inner tensor-vector products
         kr_inner = khatrirao(U[reverse(2:N)]...)
-        Rn = reshape(X, size(X, 1), :) * kr_inner
+        mul!(Rn, reshape(X, size(X, 1), :), kr_inner)
     elseif n == N
         # Just outer tensor-vector products
         kr_outer = khatrirao(U[reverse(1:N-1)]...)
-        Rn = transpose(reshape(X, prod(size(X)[1:N-1]), size(X)[N])) * kr_outer
+        mul!(Rn, transpose(reshape(X, prod(size(X)[1:N-1]), size(X)[N])), kr_outer)
     else
         for j in 1:r
             # Inner tensor-vector products

From 020a758255a7b9833b36695b165c4ddd3a433de9 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Mon, 8 Jan 2024 12:23:41 -0500
Subject: [PATCH 14/40] Try alternative Khatri-Rao def from #34

---
 src/gcp-opt.jl | 19 +++++++++++++------
 1 file changed, 13 insertions(+), 6 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 37b2daa..8233c62 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -224,14 +224,21 @@ end
 """
 function khatrirao(A::Vararg{T,N}) where {T<:Matrix,N}
     r = (only∘unique)(size.(A,2))
-    K = similar(A[1], prod(size.(A,1)), r)
-    for j in 1:r
+    #K = similar(A[1], prod(size.(A,1)), r)
+    R = ntuple(Val(N)) do k
+        dims = (ntuple(i->1,Val(N-k))..., :, ntuple(i->1,Val(k-1))..., r)
+        return reshape(A[k],dims)
+    end
+    return reshape(broadcast(*, R...),:,r)
+
+
+    #for j in 1:r
         # Quick way to avoid making two copies of `K`:
         # 1. use `reduce` to compute all but the last `kron`
         # 2. use `kron!` for the last one
         # Can/should optimize later: https://github.com/dahong67/GCPDecompositions.jl/issues/34
-        temp = reduce(kron, [view(A[i], :, j) for i in 1:N-1])
-        kron!(view(K, :, j), temp, view(A[N], :, j))
-    end
-    return K
+    #    temp = reduce(kron, [view(A[i], :, j) for i in 1:N-1])
+    #    kron!(view(K, :, j), temp, view(A[N], :, j))
+    #end
+    #return K
 end

From c31bd82cb0029eeee0a1aff8eeb4644098b54228 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Mon, 8 Jan 2024 12:59:42 -0500
Subject: [PATCH 15/40] Add benchmark suite for MTTKRP for order-4 tensors

---
 benchmark/suites/mttkrp_large.jl | 41 ++++++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)
 create mode 100644 benchmark/suites/mttkrp_large.jl

diff --git a/benchmark/suites/mttkrp_large.jl b/benchmark/suites/mttkrp_large.jl
new file mode 100644
index 0000000..4dff8b3
--- /dev/null
+++ b/benchmark/suites/mttkrp_large.jl
@@ -0,0 +1,41 @@
+module BenchmarkMTTKRP
+
+using BenchmarkTools, GCPDecompositions
+using Random
+
+const SUITE = BenchmarkGroup()
+
+# Collect setups
+const SETUPS = []
+
+## Balanced order-4 tensors
+append!(
+    SETUPS,
+    [
+        (; size = sz, rank = r, mode = n) for
+        sz in [ntuple(n -> In, 4) for In in 30:30:120], r in 30:30:180, n in 1:4
+    ],
+)
+
+## Imbalanced order-4 tensors
+append!(
+    SETUPS,
+    [
+        (; size = sz, rank = r, mode = n) for sz in [(30, 60, 100, 1000), (1000, 100, 60, 30)],
+        r in 100:100:300, n in 1:4
+    ]
+)
+
+# Generate random benchmarks
+for SETUP in SETUPS
+    Random.seed!(0)
+    X = randn(SETUP.size)
+    U = [randn(In, SETUP.rank) for In in SETUP.size]
+    SUITE["size=$(SETUP.size), rank=$(SETUP.rank), mode=$(SETUP.mode)"] = @benchmarkable(
+        GCPDecompositions.mttkrp($X, $U, $(SETUP.mode)),
+        seconds = 2,
+        samples = 5,
+    )
+end
+
+end

From 23a00812a7192fbbb45f15ed62ac539d5aa428fd Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Mon, 8 Jan 2024 13:01:54 -0500
Subject: [PATCH 16/40] Add new suite to benchmarks.jl

---
 benchmark/benchmarks.jl          | 3 ++-
 benchmark/suites/mttkrp_large.jl | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 878552b..8b433fb 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -1,7 +1,8 @@
 using BenchmarkTools
 
 # Benchmark suite modules
-const SUITE_MODULES = Dict("gcp" => :BenchmarkGCP, "mttkrp" => :BenchmarkMTTKRP)
+const SUITE_MODULES = Dict("gcp" => :BenchmarkGCP, "mttkrp" => :BenchmarkMTTKRP, 
+                            "mttkrp-large" => :BenchmarkmarkMTTKRPLarge)
 
 # Create top-level suite including only sub-suites
 # specified by ENV variable "GCP_BENCHMARK_SUITES"
diff --git a/benchmark/suites/mttkrp_large.jl b/benchmark/suites/mttkrp_large.jl
index 4dff8b3..92ff9bd 100644
--- a/benchmark/suites/mttkrp_large.jl
+++ b/benchmark/suites/mttkrp_large.jl
@@ -1,4 +1,4 @@
-module BenchmarkMTTKRP
+module BenchmarkMTTKRPLarge
 
 using BenchmarkTools, GCPDecompositions
 using Random

From c8c89d54ca168a642a9b9c50ff7af2bab013c10a Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Mon, 8 Jan 2024 13:05:55 -0500
Subject: [PATCH 17/40] Previous khatri-rao commit for comparison on new suite

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

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 8233c62..9d254bb 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -224,21 +224,9 @@ end
 """
 function khatrirao(A::Vararg{T,N}) where {T<:Matrix,N}
     r = (only∘unique)(size.(A,2))
-    #K = similar(A[1], prod(size.(A,1)), r)
-    R = ntuple(Val(N)) do k
-        dims = (ntuple(i->1,Val(N-k))..., :, ntuple(i->1,Val(k-1))..., r)
-        return reshape(A[k],dims)
+    K = similar(A[1], prod(size.(A,1)), r)
+    for j in 1:r
+        K[:, j] = reduce(kron, [view(A[i], :, j) for i in 1:N])
     end
-    return reshape(broadcast(*, R...),:,r)
-
-
-    #for j in 1:r
-        # Quick way to avoid making two copies of `K`:
-        # 1. use `reduce` to compute all but the last `kron`
-        # 2. use `kron!` for the last one
-        # Can/should optimize later: https://github.com/dahong67/GCPDecompositions.jl/issues/34
-    #    temp = reduce(kron, [view(A[i], :, j) for i in 1:N-1])
-    #    kron!(view(K, :, j), temp, view(A[N], :, j))
-    #end
-    #return K
+    return K
 end

From dc4930c24b6713db8b10868d511ff0cd9a5c2184 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Mon, 8 Jan 2024 13:21:08 -0500
Subject: [PATCH 18/40] Fix suite names

---
 benchmark/benchmarks.jl                               | 2 +-
 benchmark/suites/{mttkrp_large.jl => mttkrp-large.jl} | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)
 rename benchmark/suites/{mttkrp_large.jl => mttkrp-large.jl} (80%)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 8b433fb..4af2e86 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -2,7 +2,7 @@ using BenchmarkTools
 
 # Benchmark suite modules
 const SUITE_MODULES = Dict("gcp" => :BenchmarkGCP, "mttkrp" => :BenchmarkMTTKRP, 
-                            "mttkrp-large" => :BenchmarkmarkMTTKRPLarge)
+                            "mttkrp-large" => :BenchmarkMTTKRPLarge)
 
 # Create top-level suite including only sub-suites
 # specified by ENV variable "GCP_BENCHMARK_SUITES"
diff --git a/benchmark/suites/mttkrp_large.jl b/benchmark/suites/mttkrp-large.jl
similarity index 80%
rename from benchmark/suites/mttkrp_large.jl
rename to benchmark/suites/mttkrp-large.jl
index 92ff9bd..8b92162 100644
--- a/benchmark/suites/mttkrp_large.jl
+++ b/benchmark/suites/mttkrp-large.jl
@@ -13,7 +13,7 @@ append!(
     SETUPS,
     [
         (; size = sz, rank = r, mode = n) for
-        sz in [ntuple(n -> In, 4) for In in 30:30:120], r in 30:30:180, n in 1:4
+        sz in [ntuple(n -> In, 4) for In in 20:20:80], r in 20:20:120, n in 1:4
     ],
 )
 
@@ -21,7 +21,7 @@ append!(
 append!(
     SETUPS,
     [
-        (; size = sz, rank = r, mode = n) for sz in [(30, 60, 100, 1000), (1000, 100, 60, 30)],
+        (; size = sz, rank = r, mode = n) for sz in [(20, 40, 80, 500), (500, 80, 40, 20)],
         r in 100:100:300, n in 1:4
     ]
 )

From 0aee58d4535799c78ff4a5e01122f43c7d8b00e0 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Mon, 8 Jan 2024 13:27:02 -0500
Subject: [PATCH 19/40] Rehsape-broadcast version of khatri-rao for benchmark
 comparison

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

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 9d254bb..0c7c40a 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -222,11 +222,12 @@ end
     Computes the Khatri-Rao product (i.e., column-wise Kronecker product)
     of the matrices `A1`, `A2`, etc.
 """
-function khatrirao(A::Vararg{T,N}) where {T<:Matrix,N}
-    r = (only∘unique)(size.(A,2))
-    K = similar(A[1], prod(size.(A,1)), r)
-    for j in 1:r
-        K[:, j] = reduce(kron, [view(A[i], :, j) for i in 1:N])
+function khatrirao(A::Vararg{T,N}) where {T<:AbstractMatrix,N}
+    r = size(A[1],2)
+    # @boundscheck all(==(r),size.(A,2)) || throw(DimensionMismatch())
+    R = ntuple(Val(N)) do k
+        dims = (ntuple(i->1,Val(N-k))..., :, ntuple(i->1,Val(k-1))..., r)
+        return reshape(A[k],dims)
     end
-    return K
+    return reshape(broadcast(*, R...),:,r)
 end

From 24fda38bfbe43ccbbbae2b5bb6a1ff76023aaef6 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Tue, 9 Jan 2024 09:50:18 -0500
Subject: [PATCH 20/40] Try pre-forming inner khatri-rao for mode n != 1 and !=
 N

---
 src/gcp-opt.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 0c7c40a..10bcd1f 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -202,10 +202,10 @@ function mttkrp(X, U, n)
         kr_outer = khatrirao(U[reverse(1:N-1)]...)
         mul!(Rn, transpose(reshape(X, prod(size(X)[1:N-1]), size(X)[N])), kr_outer)
     else
+        kr_inner = khatrirao(U[reverse(n+1:N)]...)
         for j in 1:r
             # Inner tensor-vector products
-            kr_inner = reduce(kron, [view(U[i], :, j) for i in reverse(n+1:N)])
-            inner = reshape(reshape(X, Jn, Kn) * kr_inner, size(X)[1:n]) 
+            inner = reshape(reshape(X, Jn, Kn) * kr_inner[:, j], size(X)[1:n]) 
             # Outer tensor-vector products
             Jn_inner = prod(size(inner)[1:n-1])
             Kn_inner = prod(size(inner)[n:end])

From 3cf199479857bb348e37d166a041d1ba1f053cc3 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Tue, 9 Jan 2024 10:06:20 -0500
Subject: [PATCH 21/40] Pre-form outer khatri-rao for n != 1 and != N

---
 src/gcp-opt.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 10bcd1f..d62a62a 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -202,15 +202,15 @@ function mttkrp(X, U, n)
         kr_outer = khatrirao(U[reverse(1:N-1)]...)
         mul!(Rn, transpose(reshape(X, prod(size(X)[1:N-1]), size(X)[N])), kr_outer)
     else
-        kr_inner = khatrirao(U[reverse(n+1:N)]...)
+        kr_inner = khatrirao(U[reverse(n+1:N)]...) 
+        kr_outer = khatrirao(U[reverse(1:n-1)]...)
         for j in 1:r
             # Inner tensor-vector products
             inner = reshape(reshape(X, Jn, Kn) * kr_inner[:, j], size(X)[1:n]) 
             # Outer tensor-vector products
             Jn_inner = prod(size(inner)[1:n-1])
             Kn_inner = prod(size(inner)[n:end])
-            kr_outer = reduce(kron, [view(U[i], :, j) for i in reverse(1:n-1)])
-            Rn[:, j] = transpose(reshape(inner, Jn_inner, Kn_inner)) * kr_outer
+            Rn[:, j] = transpose(reshape(inner, Jn_inner, Kn_inner)) * kr_outer[:, j]
         end
     end
     return Rn

From 18e5213673e8c4481c8de1a1c13abdabdf5e8813 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Tue, 9 Jan 2024 10:28:32 -0500
Subject: [PATCH 22/40] Move inner multiply out of loop

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

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index d62a62a..c1a6f1e 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -204,13 +204,14 @@ function mttkrp(X, U, n)
     else
         kr_inner = khatrirao(U[reverse(n+1:N)]...) 
         kr_outer = khatrirao(U[reverse(1:n-1)]...)
+        inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r)) 
         for j in 1:r
             # Inner tensor-vector products
-            inner = reshape(reshape(X, Jn, Kn) * kr_inner[:, j], size(X)[1:n]) 
+            #inner = reshape(reshape(X, Jn, Kn) * kr_inner[:, j], size(X)[1:n]) 
             # Outer tensor-vector products
-            Jn_inner = prod(size(inner)[1:n-1])
-            Kn_inner = prod(size(inner)[n:end])
-            Rn[:, j] = transpose(reshape(inner, Jn_inner, Kn_inner)) * kr_outer[:, j]
+            Jn_inner = prod(size(selectdim(inner, ndims(inner), 1))[1:n-1])
+            Kn_inner = prod(size(selectdim(inner, ndims(inner), 1))[n:end])
+            Rn[:, j] = transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j]
         end
     end
     return Rn

From 2a65e1a4cc605594ea05a9db159e27c735a5ef9f Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Tue, 9 Jan 2024 10:43:37 -0500
Subject: [PATCH 23/40] Refactor Jn_inner and Kn_inner computation

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

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index c1a6f1e..62b6c6e 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -205,12 +205,9 @@ function mttkrp(X, U, n)
         kr_inner = khatrirao(U[reverse(n+1:N)]...) 
         kr_outer = khatrirao(U[reverse(1:n-1)]...)
         inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r)) 
+        Jn_inner = prod(size(inner)[1:n-1])
+        Kn_inner = prod(size(inner)[n:end-1])
         for j in 1:r
-            # Inner tensor-vector products
-            #inner = reshape(reshape(X, Jn, Kn) * kr_inner[:, j], size(X)[1:n]) 
-            # Outer tensor-vector products
-            Jn_inner = prod(size(selectdim(inner, ndims(inner), 1))[1:n-1])
-            Kn_inner = prod(size(selectdim(inner, ndims(inner), 1))[n:end])
             Rn[:, j] = transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j]
         end
     end

From b96eb46f981eee0d3cfa9880369240d977ce32cd Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Thu, 11 Jan 2024 09:39:37 -0500
Subject: [PATCH 24/40] Try .= for final result for n != 1,N case

---
 src/gcp-opt.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 62b6c6e..f95a45d 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -208,7 +208,7 @@ function mttkrp(X, U, n)
         Jn_inner = prod(size(inner)[1:n-1])
         Kn_inner = prod(size(inner)[n:end-1])
         for j in 1:r
-            Rn[:, j] = transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j]
+            Rn[:, j] .= transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j]
         end
     end
     return Rn

From b56d5d611a9c81de6bc97fe55c15d217896ff557 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Thu, 11 Jan 2024 10:19:14 -0500
Subject: [PATCH 25/40] Use reduce and hcat instead of loop

---
 src/gcp-opt.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index f95a45d..51850ac 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -207,9 +207,7 @@ function mttkrp(X, U, n)
         inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r)) 
         Jn_inner = prod(size(inner)[1:n-1])
         Kn_inner = prod(size(inner)[n:end-1])
-        for j in 1:r
-            Rn[:, j] .= transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j]
-        end
+        Rn = reduce(hcat, [transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j] for j in 1:r])
     end
     return Rn
 end

From f57a4dd77f51f6246b8ceed8e8999e0526332346 Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Thu, 11 Jan 2024 11:13:18 -0500
Subject: [PATCH 26/40] Try reusing inner

---
 src/gcp-opt.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 51850ac..b9e83c4 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -186,7 +186,6 @@ function mttkrp(X, U, n)
     (N == ndims(X) && I == size(X)) || throw(DimensionMismatch("`X` and `U` do not have matching dimensions"))
 
     # See section III-B from "Fast Alternating LS Algorithms for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al.
-    Rn = similar(U[n])
     Jn = prod(size(X)[1:n])
     Kn = prod(size(X)[n+1:end])
     # Compute tensor-vector products right to left (equations 15, 17) for each rank
@@ -195,10 +194,12 @@ function mttkrp(X, U, n)
     # n = N has no inner tensor-vector products
     if n == 1
         # Just inner tensor-vector products
+        Rn = similar(U[n])
         kr_inner = khatrirao(U[reverse(2:N)]...)
         mul!(Rn, reshape(X, size(X, 1), :), kr_inner)
     elseif n == N
         # Just outer tensor-vector products
+        Rn = similar(U[n])
         kr_outer = khatrirao(U[reverse(1:N-1)]...)
         mul!(Rn, transpose(reshape(X, prod(size(X)[1:N-1]), size(X)[N])), kr_outer)
     else
@@ -207,7 +208,8 @@ function mttkrp(X, U, n)
         inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r)) 
         Jn_inner = prod(size(inner)[1:n-1])
         Kn_inner = prod(size(inner)[n:end-1])
-        Rn = reduce(hcat, [transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j] for j in 1:r])
+        inner = reduce(hcat, [transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j] for j in 1:r])
+        return inner
     end
     return Rn
 end

From 6d33f3b8faea6b4517c518b14baa65d64752ed8a Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Thu, 11 Jan 2024 13:49:02 -0500
Subject: [PATCH 27/40] Clean up

---
 src/gcp-opt.jl | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index b9e83c4..51850ac 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -186,6 +186,7 @@ function mttkrp(X, U, n)
     (N == ndims(X) && I == size(X)) || throw(DimensionMismatch("`X` and `U` do not have matching dimensions"))
 
     # See section III-B from "Fast Alternating LS Algorithms for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al.
+    Rn = similar(U[n])
     Jn = prod(size(X)[1:n])
     Kn = prod(size(X)[n+1:end])
     # Compute tensor-vector products right to left (equations 15, 17) for each rank
@@ -194,12 +195,10 @@ function mttkrp(X, U, n)
     # n = N has no inner tensor-vector products
     if n == 1
         # Just inner tensor-vector products
-        Rn = similar(U[n])
         kr_inner = khatrirao(U[reverse(2:N)]...)
         mul!(Rn, reshape(X, size(X, 1), :), kr_inner)
     elseif n == N
         # Just outer tensor-vector products
-        Rn = similar(U[n])
         kr_outer = khatrirao(U[reverse(1:N-1)]...)
         mul!(Rn, transpose(reshape(X, prod(size(X)[1:N-1]), size(X)[N])), kr_outer)
     else
@@ -208,8 +207,7 @@ function mttkrp(X, U, n)
         inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r)) 
         Jn_inner = prod(size(inner)[1:n-1])
         Kn_inner = prod(size(inner)[n:end-1])
-        inner = reduce(hcat, [transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j] for j in 1:r])
-        return inner
+        Rn = reduce(hcat, [transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j] for j in 1:r])
     end
     return Rn
 end

From 48d965e5a2d332bff70ba7c2cefef6b933618a5a Mon Sep 17 00:00:00 2001
From: alexmul1114 <96952477+alexmul1114@users.noreply.github.com>
Date: Tue, 16 Jan 2024 13:42:58 -0500
Subject: [PATCH 28/40] Split out N=1 in Khatri-rao

---
 src/gcp-opt.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 51850ac..52a55f7 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -219,6 +219,9 @@ end
     of the matrices `A1`, `A2`, etc.
 """
 function khatrirao(A::Vararg{T,N}) where {T<:AbstractMatrix,N}
+    if N == 1
+        return A[1]
+    end
     r = size(A[1],2)
     # @boundscheck all(==(r),size.(A,2)) || throw(DimensionMismatch())
     R = ntuple(Val(N)) do k

From 4e1aee58b635928e38ff34a1f106a58a26b9ecd8 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 12:31:02 -0500
Subject: [PATCH 29/40] Run formatter

---
 benchmark/benchmarks.jl          |  7 +++++--
 benchmark/suites/mttkrp-large.jl |  6 +++---
 src/gcp-opt.jl                   | 26 ++++++++++++++++----------
 3 files changed, 24 insertions(+), 15 deletions(-)

diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 4af2e86..8ba60d3 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -1,8 +1,11 @@
 using BenchmarkTools
 
 # Benchmark suite modules
-const SUITE_MODULES = Dict("gcp" => :BenchmarkGCP, "mttkrp" => :BenchmarkMTTKRP, 
-                            "mttkrp-large" => :BenchmarkMTTKRPLarge)
+const SUITE_MODULES = Dict(
+    "gcp" => :BenchmarkGCP,
+    "mttkrp" => :BenchmarkMTTKRP,
+    "mttkrp-large" => :BenchmarkMTTKRPLarge,
+)
 
 # Create top-level suite including only sub-suites
 # specified by ENV variable "GCP_BENCHMARK_SUITES"
diff --git a/benchmark/suites/mttkrp-large.jl b/benchmark/suites/mttkrp-large.jl
index 8b92162..cce0839 100644
--- a/benchmark/suites/mttkrp-large.jl
+++ b/benchmark/suites/mttkrp-large.jl
@@ -12,8 +12,8 @@ const SETUPS = []
 append!(
     SETUPS,
     [
-        (; size = sz, rank = r, mode = n) for
-        sz in [ntuple(n -> In, 4) for In in 20:20:80], r in 20:20:120, n in 1:4
+        (; size = sz, rank = r, mode = n) for sz in [ntuple(n -> In, 4) for In in 20:20:80],
+        r in 20:20:120, n in 1:4
     ],
 )
 
@@ -23,7 +23,7 @@ append!(
     [
         (; size = sz, rank = r, mode = n) for sz in [(20, 40, 80, 500), (500, 80, 40, 20)],
         r in 100:100:300, n in 1:4
-    ]
+    ],
 )
 
 # Generate random benchmarks
diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 52a55f7..c766901 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -92,7 +92,7 @@ function _gcp(
         M0.U[k] .*= (Xnorm / M0norm)^(1 / N)
     end
     u0 = vcat(vec.(M0.U)...)
-    
+
     # Setup vectorized objective function and gradient
     vec_cutoffs = (0, cumsum(r .* size(X))...)
     vec_ranges = ntuple(k -> vec_cutoffs[k]+1:vec_cutoffs[k+1], Val(N))
@@ -182,8 +182,9 @@ end
 function mttkrp(X, U, n)
 
     # Dimensions
-    N, I, r = length(U), Tuple(size.(U, 1)), (only∘unique)(size.(U, 2))
-    (N == ndims(X) && I == size(X)) || throw(DimensionMismatch("`X` and `U` do not have matching dimensions"))
+    N, I, r = length(U), Tuple(size.(U, 1)), (only ∘ unique)(size.(U, 2))
+    (N == ndims(X) && I == size(X)) ||
+        throw(DimensionMismatch("`X` and `U` do not have matching dimensions"))
 
     # See section III-B from "Fast Alternating LS Algorithms for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al.
     Rn = similar(U[n])
@@ -202,12 +203,17 @@ function mttkrp(X, U, n)
         kr_outer = khatrirao(U[reverse(1:N-1)]...)
         mul!(Rn, transpose(reshape(X, prod(size(X)[1:N-1]), size(X)[N])), kr_outer)
     else
-        kr_inner = khatrirao(U[reverse(n+1:N)]...) 
+        kr_inner = khatrirao(U[reverse(n+1:N)]...)
         kr_outer = khatrirao(U[reverse(1:n-1)]...)
-        inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r)) 
+        inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r))
         Jn_inner = prod(size(inner)[1:n-1])
         Kn_inner = prod(size(inner)[n:end-1])
-        Rn = reduce(hcat, [transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j] for j in 1:r])
+        Rn = reduce(
+            hcat,
+            [
+                transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j] for j in 1:r
+            ],
+        )
     end
     return Rn
 end
@@ -222,11 +228,11 @@ function khatrirao(A::Vararg{T,N}) where {T<:AbstractMatrix,N}
     if N == 1
         return A[1]
     end
-    r = size(A[1],2)
+    r = size(A[1], 2)
     # @boundscheck all(==(r),size.(A,2)) || throw(DimensionMismatch())
     R = ntuple(Val(N)) do k
-        dims = (ntuple(i->1,Val(N-k))..., :, ntuple(i->1,Val(k-1))..., r)
-        return reshape(A[k],dims)
+        dims = (ntuple(i -> 1, Val(N - k))..., :, ntuple(i -> 1, Val(k - 1))..., r)
+        return reshape(A[k], dims)
     end
-    return reshape(broadcast(*, R...),:,r)
+    return reshape(broadcast(*, R...), :, r)
 end

From 357b77436599b1d21a52e8e3c7024a5149519356 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 12:35:33 -0500
Subject: [PATCH 30/40] Undo some formatting changes

---
 src/gcp-opt.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index c766901..3a51d93 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -65,6 +65,7 @@ function _gcp(
 ) where {TX,N}
     # T = promote_type(nonmissingtype(TX), Float64)
     T = Float64    # LBFGSB.jl seems to only support Float64
+
     # Compute lower bound from constraints
     lower = maximum(constraint.value for constraint in constraints; init = T(-Inf))
 
@@ -106,6 +107,7 @@ function _gcp(
         gcp_grad_U!(GU, CPD(ones(T, r), U), X, loss)
         return gu
     end
+
     # Run LBFGSB
     lbfgsopts = (; (pn => getproperty(algorithm, pn) for pn in propertynames(algorithm))...)
     u = lbfgsb(f, g!, u0; lb = fill(lower, length(u0)), lbfgsopts...)[2]
@@ -180,7 +182,6 @@ end
     section III-B.
 """
 function mttkrp(X, U, n)
-
     # Dimensions
     N, I, r = length(U), Tuple(size.(U, 1)), (only ∘ unique)(size.(U, 2))
     (N == ndims(X) && I == size(X)) ||

From 5262d61e61fcb706e4cc6283989cecf4d1e9a9a8 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 12:36:38 -0500
Subject: [PATCH 31/40] Undo some more formatting changes

To make the diff simpler for this PR.
---
 src/gcp-opt.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 3a51d93..348b888 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -183,9 +183,8 @@ end
 """
 function mttkrp(X, U, n)
     # Dimensions
-    N, I, r = length(U), Tuple(size.(U, 1)), (only ∘ unique)(size.(U, 2))
-    (N == ndims(X) && I == size(X)) ||
-        throw(DimensionMismatch("`X` and `U` do not have matching dimensions"))
+    N, I, r = length(U), Tuple(size.(U, 1)), (only∘unique)(size.(U, 2))
+    (N == ndims(X) && I == size(X)) || throw(DimensionMismatch("`X` and `U` do not have matching dimensions"))
 
     # See section III-B from "Fast Alternating LS Algorithms for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al.
     Rn = similar(U[n])

From 81be64f7b023b69c8ac7a4f2974780f8f17d7f27 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 16:13:06 -0500
Subject: [PATCH 32/40] Tweak/simplify n==N case

---
 src/gcp-opt.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 348b888..e4aafde 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -201,7 +201,7 @@ function mttkrp(X, U, n)
     elseif n == N
         # Just outer tensor-vector products
         kr_outer = khatrirao(U[reverse(1:N-1)]...)
-        mul!(Rn, transpose(reshape(X, prod(size(X)[1:N-1]), size(X)[N])), kr_outer)
+        mul!(Rn, transpose(reshape(X, :, size(X, N))), kr_outer)
     else
         kr_inner = khatrirao(U[reverse(n+1:N)]...)
         kr_outer = khatrirao(U[reverse(1:n-1)]...)

From b400538686f0d95bf54a000b20a0adcb456aa576 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 19:46:21 -0500
Subject: [PATCH 33/40] Do outer tensor-vector products in-place

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

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index e4aafde..e932ace 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -208,12 +208,13 @@ function mttkrp(X, U, n)
         inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r))
         Jn_inner = prod(size(inner)[1:n-1])
         Kn_inner = prod(size(inner)[n:end-1])
-        Rn = reduce(
-            hcat,
-            [
-                transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)) * kr_outer[:, j] for j in 1:r
-            ],
-        )
+        for j in 1:r
+            mul!(
+                view(Rn, :, j),
+                transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)),
+                view(kr_outer, :, j),
+            )
+        end
     end
     return Rn
 end

From 58b352fccda7d8a4923bdec90c4be4dbd6ddfe51 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 20:04:27 -0500
Subject: [PATCH 34/40] Add left-to-right path

---
 src/gcp-opt.jl | 40 ++++++++++++++++++++++++++++++----------
 1 file changed, 30 insertions(+), 10 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index e932ace..a5493ab 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -203,17 +203,37 @@ function mttkrp(X, U, n)
         kr_outer = khatrirao(U[reverse(1:N-1)]...)
         mul!(Rn, transpose(reshape(X, :, size(X, N))), kr_outer)
     else
-        kr_inner = khatrirao(U[reverse(n+1:N)]...)
-        kr_outer = khatrirao(U[reverse(1:n-1)]...)
-        inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r))
-        Jn_inner = prod(size(inner)[1:n-1])
-        Kn_inner = prod(size(inner)[n:end-1])
-        for j in 1:r
-            mul!(
-                view(Rn, :, j),
-                transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)),
-                view(kr_outer, :, j),
+        if prod(I[n:N]) < prod(I[1:n])  # left-to-right better
+            kr_inner = khatrirao(U[reverse(1:n-1)]...)
+            kr_outer = khatrirao(U[reverse(n+1:N)]...)
+            inner = reshape(
+                transpose(reshape(X, prod(size(X)[1:n-1]), :)) * kr_inner,
+                (size(X)[n:N]..., r),
             )
+            Jn_inner = prod(size(inner)[1:1])
+            Kn_inner = prod(size(inner)[2:end-1])
+            for j in 1:r
+                mul!(
+                    view(Rn, :, j),
+                    reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner),
+                    view(kr_outer, :, j),
+                )
+            end
+        else                            # right-to-left better
+            kr_inner = khatrirao(U[reverse(n+1:N)]...)
+            kr_outer = khatrirao(U[reverse(1:n-1)]...)
+            inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r))
+            Jn_inner = prod(size(inner)[1:n-1])
+            Kn_inner = prod(size(inner)[n:end-1])
+            for j in 1:r
+                mul!(
+                    view(Rn, :, j),
+                    transpose(
+                        reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner),
+                    ),
+                    view(kr_outer, :, j),
+                )
+            end
         end
     end
     return Rn

From 8f0ab020c631acf93bcf49687dc7f2a9b244ed33 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 21:43:30 -0500
Subject: [PATCH 35/40] Tweak high-level organization

---
 src/gcp-opt.jl | 81 +++++++++++++++++++++++++-------------------------
 1 file changed, 40 insertions(+), 41 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index a5493ab..ae11eb1 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -175,7 +175,7 @@ function _gcp(
 end
 
 """
-    mttkrp(X, U, n) -> Rn
+    mttkrp(X, U, n)
     
     Algorithm for computing one mode of MTTKRP is from "Fast Alternating LS Algorithms
     for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al., specifically
@@ -186,57 +186,56 @@ function mttkrp(X, U, n)
     N, I, r = length(U), Tuple(size.(U, 1)), (only∘unique)(size.(U, 2))
     (N == ndims(X) && I == size(X)) || throw(DimensionMismatch("`X` and `U` do not have matching dimensions"))
 
-    # See section III-B from "Fast Alternating LS Algorithms for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al.
-    Rn = similar(U[n])
-    Jn = prod(size(X)[1:n])
-    Kn = prod(size(X)[n+1:end])
-    # Compute tensor-vector products right to left (equations 15, 17) for each rank
+    # Allocate output array G
+    G = similar(U[n])
 
-    # Special cases are n = 1 and n = N (n = 1 has no outer tensor-vector products),
-    # n = N has no inner tensor-vector products
+    # Choose appropriate multiplication order:
+    # + n == 1: no splitting required
+    # + n == N: no splitting required
+    # + 1 < n < N: better to multiply "bigger" side out first
+    #   + prod(I[1:n]) > prod(I[n:N]): better to multiply left-to-right
+    #   + prod(I[1:n]) < prod(I[n:N]): better to multiply right-to-left
     if n == 1
         # Just inner tensor-vector products
         kr_inner = khatrirao(U[reverse(2:N)]...)
-        mul!(Rn, reshape(X, size(X, 1), :), kr_inner)
+        mul!(G, reshape(X, size(X, 1), :), kr_inner)
     elseif n == N
         # Just outer tensor-vector products
         kr_outer = khatrirao(U[reverse(1:N-1)]...)
-        mul!(Rn, transpose(reshape(X, :, size(X, N))), kr_outer)
+        mul!(G, transpose(reshape(X, :, size(X, N))), kr_outer)
+    elseif prod(I[1:n]) > prod(I[n:N])
+        kr_inner = khatrirao(U[reverse(1:n-1)]...)
+        kr_outer = khatrirao(U[reverse(n+1:N)]...)
+        inner = reshape(
+            transpose(reshape(X, prod(size(X)[1:n-1]), :)) * kr_inner,
+            (size(X)[n:N]..., r),
+        )
+        Jn_inner = prod(size(inner)[1:1])
+        Kn_inner = prod(size(inner)[2:end-1])
+        for j in 1:r
+            mul!(
+                view(G, :, j),
+                reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner),
+                view(kr_outer, :, j),
+            )
+        end
     else
-        if prod(I[n:N]) < prod(I[1:n])  # left-to-right better
-            kr_inner = khatrirao(U[reverse(1:n-1)]...)
-            kr_outer = khatrirao(U[reverse(n+1:N)]...)
-            inner = reshape(
-                transpose(reshape(X, prod(size(X)[1:n-1]), :)) * kr_inner,
-                (size(X)[n:N]..., r),
+        kr_inner = khatrirao(U[reverse(n+1:N)]...)
+        kr_outer = khatrirao(U[reverse(1:n-1)]...)
+        Jn = prod(size(X)[1:n])
+        Kn = prod(size(X)[n+1:end])
+        inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r))
+        Jn_inner = prod(size(inner)[1:n-1])
+        Kn_inner = prod(size(inner)[n:end-1])
+        for j in 1:r
+            mul!(
+                view(G, :, j),
+                transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)),
+                view(kr_outer, :, j),
             )
-            Jn_inner = prod(size(inner)[1:1])
-            Kn_inner = prod(size(inner)[2:end-1])
-            for j in 1:r
-                mul!(
-                    view(Rn, :, j),
-                    reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner),
-                    view(kr_outer, :, j),
-                )
-            end
-        else                            # right-to-left better
-            kr_inner = khatrirao(U[reverse(n+1:N)]...)
-            kr_outer = khatrirao(U[reverse(1:n-1)]...)
-            inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r))
-            Jn_inner = prod(size(inner)[1:n-1])
-            Kn_inner = prod(size(inner)[n:end-1])
-            for j in 1:r
-                mul!(
-                    view(Rn, :, j),
-                    transpose(
-                        reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner),
-                    ),
-                    view(kr_outer, :, j),
-                )
-            end
         end
     end
-    return Rn
+    return G
 end
 
 """

From 9ade08cc1737271e01acb81f2af4e7fd30a82293 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 22:14:36 -0500
Subject: [PATCH 36/40] Tidy up / simplify implementations

---
 src/gcp-opt.jl | 43 ++++++++++++++++++-------------------------
 1 file changed, 18 insertions(+), 25 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index ae11eb1..c66ccb2 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -196,42 +196,35 @@ function mttkrp(X, U, n)
     #   + prod(I[1:n]) > prod(I[n:N]): better to multiply left-to-right
     #   + prod(I[1:n]) < prod(I[n:N]): better to multiply right-to-left
     if n == 1
-        # Just inner tensor-vector products
-        kr_inner = khatrirao(U[reverse(2:N)]...)
-        mul!(G, reshape(X, size(X, 1), :), kr_inner)
+        mul!(G, reshape(X, I[1], :), khatrirao(U[reverse(2:N)]...))
     elseif n == N
-        # Just outer tensor-vector products
-        kr_outer = khatrirao(U[reverse(1:N-1)]...)
-        mul!(G, transpose(reshape(X, :, size(X, N))), kr_outer)
+        mul!(G, transpose(reshape(X, :, I[N])), khatrirao(U[reverse(1:N-1)]...))
     elseif prod(I[1:n]) > prod(I[n:N])
-        kr_inner = khatrirao(U[reverse(1:n-1)]...)
-        kr_outer = khatrirao(U[reverse(n+1:N)]...)
-        inner = reshape(
-            transpose(reshape(X, prod(size(X)[1:n-1]), :)) * kr_inner,
-            (size(X)[n:N]..., r),
-        )
-        Jn_inner = prod(size(inner)[1:1])
-        Kn_inner = prod(size(inner)[2:end-1])
+        # Inner multiplication: left side
+        kr_left = khatrirao(U[reverse(1:n-1)]...)
+        L = reshape(transpose(reshape(X, :, prod(I[n:N]))) * kr_left, (I[n:N]..., r))
+
+        # Outer multiplication: right side
+        kr_right = khatrirao(U[reverse(n+1:N)]...)
         for j in 1:r
             mul!(
                 view(G, :, j),
-                reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner),
-                view(kr_outer, :, j),
+                reshape(selectdim(L, ndims(L), j), I[n], :),
+                view(kr_right, :, j),
             )
         end
     else
-        kr_inner = khatrirao(U[reverse(n+1:N)]...)
-        kr_outer = khatrirao(U[reverse(1:n-1)]...)
-        Jn = prod(size(X)[1:n])
-        Kn = prod(size(X)[n+1:end])
-        inner = reshape(reshape(X, Jn, Kn) * kr_inner, (size(X)[1:n]..., r))
-        Jn_inner = prod(size(inner)[1:n-1])
-        Kn_inner = prod(size(inner)[n:end-1])
+        # Inner multiplication: right side
+        kr_right = khatrirao(U[reverse(n+1:N)]...)
+        R = reshape(reshape(X, prod(I[1:n]), :) * kr_right, (I[1:n]..., r))
+
+        # Outer multiplication: left side
+        kr_left = khatrirao(U[reverse(1:n-1)]...)
         for j in 1:r
             mul!(
                 view(G, :, j),
-                transpose(reshape(selectdim(inner, ndims(inner), j), Jn_inner, Kn_inner)),
-                view(kr_outer, :, j),
+                transpose(reshape(selectdim(R, ndims(R), j), :, I[n])),
+                view(kr_left, :, j),
             )
         end
     end

From 5322a3af9cbbc85cfd103d9bd095df399c017b2d Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 22:28:28 -0500
Subject: [PATCH 37/40] Uncomment dimension check in khatrirao

---
 src/gcp-opt.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index c66ccb2..b26aad1 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -238,11 +238,14 @@ end
     of the matrices `A1`, `A2`, etc.
 """
 function khatrirao(A::Vararg{T,N}) where {T<:AbstractMatrix,N}
+    # Special case: N = 1
     if N == 1
         return A[1]
     end
+
+    # General case: N > 1
     r = size(A[1], 2)
-    # @boundscheck all(==(r),size.(A,2)) || throw(DimensionMismatch())
+    all(==(r),size.(A,2)) || throw(DimensionMismatch())
     R = ntuple(Val(N)) do k
         dims = (ntuple(i -> 1, Val(N - k))..., :, ntuple(i -> 1, Val(k - 1))..., r)
         return reshape(A[k], dims)

From 1db9d0f31098e8e69a11d7270f0e0b32a2fc6011 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 22:39:56 -0500
Subject: [PATCH 38/40] Edit docstrings

---
 src/gcp-opt.jl | 21 +++++++++++++--------
 1 file changed, 13 insertions(+), 8 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index b26aad1..8343d6e 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -175,11 +175,16 @@ function _gcp(
 end
 
 """
-    mttkrp(X, U, n)
-    
-    Algorithm for computing one mode of MTTKRP is from "Fast Alternating LS Algorithms
-    for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al., specifically
-    section III-B.
+    mttkrp(X, (U1, U2, ..., UN), n)
+
+Compute the Matricized Tensor Times Khatri-Rao Product (MTTKRP)
+of an N-way tensor X with the matrices U1, U2, ..., UN along mode n.
+
+Algorithm is based on Section III-B of the paper:
+> **Fast Alternating LS Algorithms for High Order CANDECOMP/PARAFAC Tensor Factorizations**.
+> Anh-Huy Phan, Petr Tichavský, Andrzej Cichocki.
+> *IEEE Transactions on Signal Processing*, 2013.
+> DOI: 10.1109/TSP.2013.2269903
 """
 function mttkrp(X, U, n)
     # Dimensions
@@ -233,9 +238,9 @@ end
 
 """
     khatrirao(A1, A2, ...)
-    
-    Computes the Khatri-Rao product (i.e., column-wise Kronecker product)
-    of the matrices `A1`, `A2`, etc.
+
+Compute the Khatri-Rao product (i.e., the column-wise Kronecker product)
+of the matrices `A1`, `A2`, etc.
 """
 function khatrirao(A::Vararg{T,N}) where {T<:AbstractMatrix,N}
     # Special case: N = 1

From 684af2c3bcb54a34a15bb4b52bfea82a2f76b353 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 22:45:56 -0500
Subject: [PATCH 39/40] Reverse the order of the sizes to test the
 left-to-right MTTKRP path

---
 test/items/gcp-opt.jl | 26 +++++++++++++-------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/test/items/gcp-opt.jl b/test/items/gcp-opt.jl
index 7dba4e4..4b5086c 100644
--- a/test/items/gcp-opt.jl
+++ b/test/items/gcp-opt.jl
@@ -35,7 +35,7 @@ end
 @testitem "LeastSquaresLoss" begin
     using Random
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(ones(r), rand.(sz, r))
         X = [M[I] for I in CartesianIndices(size(M))]
@@ -55,7 +55,7 @@ end
 @testitem "NonnegativeLeastSquaresLoss" begin
     using Random
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(ones(r), rand.(sz, r))
         X = [M[I] for I in CartesianIndices(size(M))]
@@ -73,7 +73,7 @@ end
     using Random
     using Distributions
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(fill(10.0, r), rand.(sz, r))
         X = [rand(Poisson(M[I])) for I in CartesianIndices(size(M))]
@@ -100,7 +100,7 @@ end
     using Random
     using Distributions
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(ones(r), randn.(sz, r))
         X = [rand(Poisson(exp(M[I]))) for I in CartesianIndices(size(M))]
@@ -127,7 +127,7 @@ end
     using Random
     using Distributions
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(ones(r), rand.(sz, r))
         k = 1.5
@@ -155,7 +155,7 @@ end
     using Random
     using Distributions
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(ones(r), rand.(sz, r))
         X = [rand(Rayleigh(M[I]/(sqrt(pi/2)))) for I in CartesianIndices(size(M))]
@@ -182,7 +182,7 @@ end
     using Random
     using Distributions
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(ones(r), rand.(sz, r))
         X = [rand(Bernoulli(M[I]/(M[I] + 1))) for I in CartesianIndices(size(M))]
@@ -209,7 +209,7 @@ end
     using Random
     using Distributions
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(ones(r), rand.(sz, r))
         X = [rand(Bernoulli(exp(M[I])/(exp(M[I]) + 1))) for I in CartesianIndices(size(M))]
@@ -236,7 +236,7 @@ end
     using Random
     using Distributions
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(ones(r), rand.(sz, r))
         num_failures = 5
@@ -265,7 +265,7 @@ end
     using Random
     using Distributions
 
-    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+    @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
         Random.seed!(0)
         M = CPD(ones(r), rand.(sz, r))
         X = [M[I] for I in CartesianIndices(size(M))]
@@ -294,7 +294,7 @@ end
     using Random
     using Distributions
 
-    @testset "size(X)=$sz, rank(X)=$r, β" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2, β in [0, 0.5, 1]
+    @testset "size(X)=$sz, rank(X)=$r, β" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2, β in [0, 0.5, 1]
         Random.seed!(0)
         M = CPD(ones(r), rand.(sz, r))
         # May want to consider other distributions depending on value of β
@@ -342,7 +342,7 @@ end
     using Random, Distributions, IntervalSets
 
     @testset "Least Squares" begin
-        @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+        @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
             Random.seed!(0)
             M = CPD(ones(r), randn.(sz, r))
             X = [M[I] for I in CartesianIndices(size(M))]
@@ -366,7 +366,7 @@ end
     end
 
     @testset "Poisson" begin
-        @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (30, 40, 50)], r in 1:2
+        @testset "size(X)=$sz, rank(X)=$r" for sz in [(15, 20, 25), (50, 40, 30)], r in 1:2
             Random.seed!(0)
             M = CPD(fill(10.0, r), rand.(sz, r))
             X = [rand(Poisson(M[I])) for I in CartesianIndices(size(M))]

From 9159ef6af1159ff58ba641339bca4298da6b7cc1 Mon Sep 17 00:00:00 2001
From: David Hong <hong@udel.edu>
Date: Thu, 29 Feb 2024 22:56:43 -0500
Subject: [PATCH 40/40] Add back simpler Khatri-Rao

To save improvements of this function for another more focused PR.
---
 src/gcp-opt.jl | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/gcp-opt.jl b/src/gcp-opt.jl
index 8343d6e..29c62cb 100644
--- a/src/gcp-opt.jl
+++ b/src/gcp-opt.jl
@@ -249,11 +249,10 @@ function khatrirao(A::Vararg{T,N}) where {T<:AbstractMatrix,N}
     end
 
     # General case: N > 1
-    r = size(A[1], 2)
-    all(==(r),size.(A,2)) || throw(DimensionMismatch())
-    R = ntuple(Val(N)) do k
-        dims = (ntuple(i -> 1, Val(N - k))..., :, ntuple(i -> 1, Val(k - 1))..., r)
-        return reshape(A[k], dims)
+    r = (only ∘ unique)(size.(A, 2))
+    K = similar(A[1], prod(size.(A, 1)), r)
+    for j in 1:r
+        K[:, j] = reduce(kron, [view(A[i], :, j) for i in 1:N])
     end
-    return reshape(broadcast(*, R...), :, r)
+    return K
 end