CPD Normalize Function

docs/src/man/main.md

@@ -9,6 +9,8 @@ GCPDecompositions
 gcp
 CPD
 ncomps
+normalizecomps
+normalizecomps!
 GCPDecompositions.default_constraints
 GCPDecompositions.default_algorithm
 GCPDecompositions.default_init

src/GCPDecompositions.jl

@@ -12,7 +12,7 @@ using Random: default_rng
 
 # Exports
 export CPD
-export ncomps
+export ncomps, normalizecomps, normalizecomps!
 export gcp
 export GCPLosses, GCPConstraints, GCPAlgorithms
 

src/cpd.jl

@@ -93,3 +93,32 @@ function norm2(M::CPD{T,N}) where {T,N}
     V = reduce(.*, M.U[i]'M.U[i] for i in 1:N)
     return sqrt(abs(M.λ' * V * M.λ))
 end
+
+"""
+    normalizecomps(M::CPD, p::Real = 2)
+
+Normalize the components of `M` so that the columns of all its factor matrices
+all have `p`-norm equal to unity, i.e., `norm(M.U[k][:, j], p) == 1` for all
+`k ∈ 1:ndims(M)` and `j ∈ 1:ncomps(M)`. The excess weight is absorbed into `M.λ`.
+
+See also: `normalizecomps!`.
+"""
+normalizecomps(M::CPD, p::Real = 2) = normalizecomps!(deepcopy(M), p)
+
+"""
+    normalizecomps!(M::CPD, p::Real = 2)
+
+Normalize the components of `M` in-place so that the columns of all its factor matrices
+all have `p`-norm equal to unity, i.e., `norm(M.U[k][:, j], p) == 1` for all
+`k ∈ 1:ndims(M)` and `j ∈ 1:ncomps(M)`. The excess weight is absorbed into `M.λ`.
+
+See also: `normalizecomps`.
+"""
+function normalizecomps!(M::CPD{T,N}, p::Real = 2) where {T,N}
+    for k in 1:N
+        norms = mapslices(Base.Fix2(norm, p), M.U[k]; dims = 1)
+        M.U[k] ./= norms
+        M.λ .*= dropdims(norms; dims = 1)
+    end
+    return M
+end

test/items/cpd.jl

@@ -157,3 +157,41 @@ end
               (sum(m -> abs(m)^3, M[I] for I in CartesianIndices(size(M))))^(1 / 3)
     end
 end
+
+@testitem "normalizecomps" begin
+    using LinearAlgebra
+
+    @testset "K=$K" for K in 1:3
+        T = Float64
+        λfull = T[1, 100, 10000]
+        U1full, U2full, U3full = T[1 2 3; 4 5 6], T[-1 2 1], T[1 2 3; 4 5 6; 7 8 9]
+        λ = λfull[1:K]
+        U1, U2, U3 = U1full[:, 1:K], U2full[:, 1:K], U3full[:, 1:K]
+
+        @testset "p=$p" for p in [1, 2, Inf]
+            M = CPD(λ, (U1, U2, U3))
+            Mback = deepcopy(M)
+            Mnorm = normalizecomps(M, p)
+
+            # Check for mutation
+            @test M.λ == Mback.λ
+            @test M.U == Mback.U
+
+            # Check factors
+            @test all(1:ndims(Mnorm)) do k
+                all(1:ncomps(Mnorm)) do j
+                    return norm(Mnorm.U[k][:, j], p) ≈ 1.0
+                end
+            end
+
+            # Check weights
+            scalings = dropdims.(mapslices.(x -> norm(x, p), M.U; dims = 1); dims = 1)
+            @test Mnorm.λ ≈ M.λ .* reduce(.*, scalings)
+
+            # Check in-place version
+            normalizecomps!(M, p)
+            @test M.λ == Mnorm.λ
+            @test M.U == Mnorm.U
+        end
+    end
+end