Merge pull request #49 from dahong67/dahong67/issue48

Project.toml

@@ -10,6 +10,7 @@ IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LossFunctions = "30fc2ffe-d236-52d8-8643-a9d8f7c094a7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [weakdeps]
 LossFunctions = "30fc2ffe-d236-52d8-8643-a9d8f7c094a7"
@@ -24,4 +25,5 @@ IntervalSets = "0.7.7"
 LBFGSB = "0.4.1"
 LinearAlgebra = "1.6"
 LossFunctions = "0.11.1"
+Random = "1.6"
 julia = "1.6"

src/GCPDecompositions.jl

@@ -8,6 +8,7 @@ import Base: ndims, size, show, summary
 import Base: getindex
 import LinearAlgebra: norm
 using IntervalSets: Interval
+using Random: default_rng
 
 # Exports
 export CPD
@@ -30,7 +31,8 @@ end
 """
     gcp(X::Array, r, loss = GCPLosses.LeastSquaresLoss();
         constraints = default_constraints(loss),
-        algorithm = default_algorithm(X, r, loss, constraints))
+        algorithm = default_algorithm(X, r, loss, constraints),
+        init = default_init(X, r, loss, constraints, algorithm))
 
 Compute an approximate rank-`r` CP decomposition of the tensor `X`
 with respect to the loss function `loss` and return a `CPD` object.
@@ -55,7 +57,8 @@ gcp(
     loss = GCPLosses.LeastSquaresLoss();
     constraints = default_constraints(loss),
     algorithm = default_algorithm(X, r, loss, constraints),
-) = GCPAlgorithms._gcp(X, r, loss, constraints, algorithm)
+    init = default_init(X, r, loss, constraints, algorithm),
+) = GCPAlgorithms._gcp(X, r, loss, constraints, algorithm, init)
 
 # Defaults
 
@@ -95,4 +98,31 @@ default_algorithm(
 ) = GCPAlgorithms.ALS()
 default_algorithm(X, r, loss, constraints) = GCPAlgorithms.LBFGSB()
 
+"""
+    default_init([rng=default_rng()], X, r, loss, constraints, algorithm)
+
+Return a default initialization for the data tensor `X`, rank `r`,
+loss function `loss`, tuple of constraints `constraints`, and
+algorithm `algorithm`, using the random number generator `rng` if needed.
+
+See also: `gcp`.
+"""
+default_init(X, r, loss, constraints, algorithm) =
+    default_init(default_rng(), X, r, loss, constraints, algorithm)
+function default_init(rng, X, r, loss, constraints, algorithm)
+    # Generate CPD with random factors
+    T, N = nonmissingtype(eltype(X)), ndims(X)
+    T = promote_type(T, Float64)
+    M = CPD(ones(T, r), rand.(rng, T, size(X), r))
+
+    # Normalize
+    Mnorm = norm(M)
+    Xnorm = sqrt(sum(abs2, skipmissing(X)))
+    for k in Base.OneTo(N)
+        M.U[k] .*= (Xnorm / Mnorm)^(1 / N)
+    end
+
+    return M
+end
+
 end

src/gcp-algorithms/als.jl

@@ -20,17 +20,10 @@ function _gcp(
     loss::GCPLosses.LeastSquaresLoss,
     constraints::Tuple{},
     algorithm::GCPAlgorithms.ALS,
+    init,
 ) where {TX<:Real,N}
-    T = promote_type(TX, Float64)
-
-    # Random initialization
-    M0 = CPD(ones(T, r), rand.(T, size(X), r))
-    M0norm = norm(M0)
-    Xnorm = sqrt(sum(abs2, skipmissing(X)))
-    for k in Base.OneTo(N)
-        M0.U[k] .*= (Xnorm / M0norm)^(1 / N)
-    end
-    M = deepcopy(M0)
+    # Initialization
+    M = deepcopy(init)
 
     # Pre-allocate MTTKRP buffers
     mttkrp_buffers = ntuple(n -> create_mttkrp_buffer(X, M.U, n), N)

src/gcp-algorithms/lbfgsb.jl

@@ -37,6 +37,7 @@ function _gcp(
     loss,
     constraints::Tuple{Vararg{GCPConstraints.LowerBound}},
     algorithm::GCPAlgorithms.LBFGSB,
+    init,
 ) where {TX,N}
     # T = promote_type(nonmissingtype(TX), Float64)
     T = Float64    # LBFGSB.jl seems to only support Float64
@@ -60,13 +61,8 @@ function _gcp(
         )
     end
 
-    # Random initialization
-    M0 = CPD(ones(T, r), rand.(T, size(X), r))
-    M0norm = sqrt(sum(abs2, M0[I] for I in CartesianIndices(size(M0))))
-    Xnorm = sqrt(sum(abs2, skipmissing(X)))
-    for k in Base.OneTo(N)
-        M0.U[k] .*= (Xnorm / M0norm)^(1 / N)
-    end
+    # Initialization
+    M0 = deepcopy(init)
     u0 = vcat(vec.(M0.U)...)
 
     # Setup vectorized objective function and gradient