clean outer_product and cpproduct

Author: njericha

src/BlockTensorFactorization.jl

@@ -12,7 +12,7 @@ export SuperDiagonal, abs_randn, all_recursive, eachfibre, getnotindex, geomean,
 export d_dx, d2_dx2, curvature, standard_curvature
 
 #include("./tensorproducts.jl")
-export nmp, nmode_product, mtt, slicewise_dot, tuckerproduct, cpproduct
+export nmp, nmode_product, mtt, slicewise_dot, tuckerproduct, cpproduct, outer_product
 export ×₁, ×₂, ×₃, ×₄, ×₅, ×₆, ×₇, ×₈, ×₉
 export ⋅₁, ⋅₂, ⋅₃, ⋅₄, ⋅₅, ⋅₆, ⋅₇, ⋅₈, ⋅₉
 

src/Core/Core.jl

@@ -22,7 +22,7 @@ include("./curvaturetools.jl")
 export d_dx, d2_dx2, curvature, standard_curvature
 
 include("./tensorproducts.jl")
-export nmp, nmode_product, mtt, slicewise_dot, tuckerproduct, cpproduct
+export nmp, nmode_product, mtt, slicewise_dot, tuckerproduct, cpproduct, outer_product
 export ×₁, ×₂, ×₃, ×₄, ×₅, ×₆, ×₇, ×₈, ×₉
 export ⋅₁, ⋅₂, ⋅₃, ⋅₄, ⋅₅, ⋅₆, ⋅₇, ⋅₈, ⋅₉
 

src/Core/decomposition.jl

@@ -441,7 +441,6 @@ end
 factors(CPD::CPDecomposition) = CPD.factors
 array(CPD::CPDecomposition) = cpproduct(factors(CPD))
 frozen(CPD::CPDecomposition) = CPD.frozen
-vector_outer(v) = reshape(kron(reverse(v)...),length.(v))
 eachfactorindex(CPD::CPDecomposition) = 1:nfactors(CPD) # unlike other AbstractTucker's, back to 1 based since there's only matrix factors
 isfrozen(CPD::CPDecomposition, n::Integer) = n == 0 ? true : frozen(CPD)[n] # similar to eachfactorindex
 

src/Core/factorize.jl

@@ -153,7 +153,7 @@ function default_kwargs(Y; kwargs...)
 
 	get!(kwargs, :rank) do # Can also be a tuple. For example, Tucker rank could be (1, 2, 3) for an order 3 array Y
 		isnothing(kwargs[:decomposition]) ? error("`rank_detect_factorize` should be called if no initial decomposition nor the rank is provided.") : rankof(kwargs[:decomposition])
-	end
+	end # TODO add a check that the keyword rank is compatible with the model and decomposition
 	get!(kwargs, :init) do
 		isnonnegative(Y) ? abs_randn : randn
 	end

src/Core/tensorproducts.jl

@@ -195,6 +195,25 @@ function _gettuckerindex(core, matrices, I)
     return Y
 end
 
+"""
+    outer_product(vectors)
+    outer_product(vectors...)
+
+Outer product of a collection of vectors.
+
+For example,
+
+`outer_product(u, v) == u * v'`
+
+and
+
+`outer_product(u, v, w)[i, j, k] == u[i] * v[j] * w[k]`.
+
+Returned array will have same number dimensions as the length of the collection.
+"""
+outer_product(vectors) = reshape(kron(reverse(vectors)...),length.(vectors))
+outer_product(vectors...) = outer_product(vectors)
+
 """
     cpproduct((A, B, C, ...))
     cpproduct(A, B, C, ...)
@@ -205,7 +224,7 @@ Example
 -------
 cpproduct(A, B, C) == @einsum T[i, j, k] := A[i, r] * B[j, r] * C[k, r]
 """
-cpproduct(matrices) = mapreduce(vector_outer, +, zip((eachcol.(matrices))...))
+cpproduct(matrices) = mapreduce(outer_product, +, zip((eachcol.(matrices))...))
 cpproduct(matrices...) = cpproduct(matrices)
 
 """

test/runtests.jl

@@ -97,6 +97,18 @@ const VERBOSE = true
         B = randn(10, 20)
         @test slicewise_dot(A, A) ≈ A*A' # test this separately since it uses a different routine when the argument is the same
         @test slicewise_dot(A, B) ≈ A*B'
+
+        # CP product
+        R = 2
+        A = randn(4, R)
+        B = randn(5, R)
+        C = randn(6, R)
+        Y = CPDecomposition((A, B, C))
+        Y_array = array(Y)
+
+        @test cpproduct(A, B, C) ≈ Y_array
+        @test Y ≈ Y_array
+        @test Y[1, 2, 3] ≈ sum(A[1, r] * B[2, r] * C[3, r] for r in 1:R)
     end
 
 @testset "Constraints" begin
@@ -628,6 +640,29 @@ end
     @testset "CPFactorization" begin
         fact = BlockTensorFactorization.factorize
 
+        # Regular run of CPDecomposition
+        A = randn(5, 2)
+        B = randn(5, 2)
+        C = randn(5, 2)
+        Y = CPDecomposition((A, B, C))
+        Y = array(Y)
+
+        decomposition, stats, kwargs = fact(Y;
+            rank=2,
+            model=CPDecomposition,
+            momentum=false,
+            #tolerance=(1, 0.045),
+            tolerance=0.045,
+            converged=RelativeError,
+            #converged=(GradientNNCone, RelativeError),
+            constrain_init=false,
+            maxiter=1,
+            #constraints=nonnegative!,
+            stats=[Iteration, ObjectiveValue, GradientNNCone, RelativeError, PrintStats]
+        );
+
+
+
         # Semi-interesting run of CPDecomposition
         N = 100
         R = 5