Various ICA fixes

src/ica.jl

@@ -1,10 +1,8 @@
 # Independent Component Analysis
 
-import Printf
-
 #### FastICA type
 
-mutable struct ICA{T<:Real}
+struct ICA{T<:Real}
     mean::Vector{T}   # mean vector, of length m (or empty to indicate zero mean)
     W::Matrix{T}      # component coefficient matrix, of size (m, k)
 end
@@ -61,8 +59,7 @@ function fastica!(W::DenseMatrix{T},      # initialized component matrix, size (
                   X::DenseMatrix{T},      # (whitened) observation sample matrix, size(m, n)
                   fun::ICAGDeriv{T},      # approximate neg-entropy functor
                   maxiter::Int,           # maximum number of iterations
-                  tol::Real,              # convergence tolerance
-                  verbose::Bool) where T<:Real          # whether to show iterative info
+                  tol::Real) where T<:Real# convergence tolerance
 
     # argument checking
     m = size(W, 1)
@@ -71,10 +68,7 @@ function fastica!(W::DenseMatrix{T},      # initialized component matrix, size (
     n = size(X, 2)
     k <= min(m, n) || throw(DimensionMismatch("k must not exceed min(m, n)."))
 
-    if verbose
-        Printf.@printf("FastICA Algorithm (m = %d, n = %d, k = %d)\n", m, n, k)
-        println("============================================")
-    end
+    @debug "FastICA Algorithm" m=m n=n k=k
 
     # pre-allocated storage
     Wp = Matrix{T}(undef, m, k)    # to store previous version of W
@@ -126,41 +120,27 @@ function fastica!(W::DenseMatrix{T},      # initialized component matrix, size (
         # symmetric decorrelation: W <- W * (W'W)^{-1/2}
         copyto!(W, W * _invsqrtm!(W'W))
 
-        # compare with Wp
-        chg = 0.0
-        for j = 1:k
-            s = 0.0
-            w = view(W,:,j)
-            wp = view(Wp,:,j)
-            for i = 1:m
-                s += abs(w[i] - wp[i])
-            end
-            if s > chg
-                chg = s
-            end
-        end
+        # compare with Wp to evaluate a conversion change
+        chg = maximum(abs.(abs.(diag(W*Wp')) .- 1.0))
         converged = (chg < tol)
 
-        if verbose
-            Printf.@printf("Iter %4d:  change = %.6e\n", t, chg)
-        end
+        @debug "Iteration $t" change=chg tolerance=tol
     end
     converged || throw(ConvergenceException(maxiter, chg, oftype(chg, tol)))
     return W
 end
 
 #### interface function
 
-function fit(::Type{ICA}, X::DenseMatrix{T},                # sample matrix, size (m, n)
+function fit(::Type{ICA}, X::AbstractMatrix{T},                # sample matrix, size (m, n)
                           k::Int;                           # number of independent components
                           alg::Symbol=:fastica,             # choice of algorithm
                           fun::ICAGDeriv=icagfun(:tanh, T), # approx neg-entropy functor
                           do_whiten::Bool=true,             # whether to perform pre-whitening
                           maxiter::Integer=100,             # maximum number of iterations
                           tol::Real=1.0e-6,                 # convergence tolerance
                           mean=nothing,                     # pre-computed mean
-                          winit::Matrix{T}=zeros(T,0,0),    # init guess of W, size (m, k)
-                          verbose::Bool=false) where T<:Real             # whether to display iterations
+                          winit::Matrix{T}=zeros(T,0,0)) where T<:Real  # init guess of W, size (m, k)
 
     # check input arguments
     m, n = size(X)
@@ -175,7 +155,7 @@ function fit(::Type{ICA}, X::DenseMatrix{T},                # sample matrix, siz
     mv = preprocess_mean(X, mean)
     Z::Matrix{T} = centralize(X, mv)
 
-    W0= zeros(T, 0,0)    # whitening matrix
+    W0= zeros(T, 0, 0)  # whitening matrix
     if do_whiten
         C = rmul!(Z * transpose(Z), 1.0 / (n - 1))
         Efac = eigen(C)
@@ -189,7 +169,7 @@ function fit(::Type{ICA}, X::DenseMatrix{T},                # sample matrix, siz
     W = (isempty(winit) ? randn(T, size(Z,1), k) : copy(winit))
 
     # invoke core algorithm
-    fastica!(W, Z, fun, maxiter, tol, verbose)
+    fastica!(W, Z, fun, maxiter, tol)
 
     # construct model
     if do_whiten

test/ica.jl

@@ -5,100 +5,116 @@ import Statistics: mean, cov
 import Random
 import StatsBase
 
-Random.seed!(15678)
-
-## icagfun
-
-f = icagfun(:tanh)
-u, v = evaluate(f, 1.5)
-@test u ≈ 0.905148253644866438242
-@test v ≈ 0.180706638923648530597
-
-f = icagfun(:tanh, Float32)
-u, v = evaluate(f, 1.5f0)
-@test u ≈ 0.90514827f0
-@test v ≈ 0.18070662f0
-
-f = icagfun(:tanh, 1.5)
-u, v = evaluate(f, 1.2)
-@test u ≈ 0.946806012846268289646
-@test v ≈ 0.155337561057228069719
-
-f = icagfun(:tanh, 1.5f0)
-u, v = evaluate(f, 1.2f0)
-@test u ≈ 0.94680610f0
-@test v ≈ 0.15533754f0
-
-f = icagfun(:gaus)
-u, v = evaluate(f, 1.5)
-@test u ≈ 0.486978701037524594696
-@test v ≈ -0.405815584197937162246
-
-f = icagfun(:gaus, Float32)
-u, v = evaluate(f, 1.5f0)
-@test u ≈ 0.4869787f0
-@test v ≈ -0.40581557f0
-
-
-## data
-
-# sources
-n = 1000
-k = 3
-m = 8
-
-t = range(0.0, step=10.0, length=n)
-s1 = sin.(t * 2)
-s2 = s2 = 1.0 .- 2.0 * Bool[isodd(floor(Int, x / 3)) for x in t]
-s3 = Float64[mod(x, 5.0) for x in t]
-
-s1 += 0.1 * randn(n)
-s2 += 0.1 * randn(n)
-s3 += 0.1 * randn(n)
-
-S = hcat(s1, s2, s3)'
-@assert size(S) == (k, n)
-
-A = randn(m, k)
-
-X = A * S
-mv = vec(mean(X, dims=2))
-@assert size(X) == (m, n)
-C = cov(X, dims=2)
-
-# FastICA
-
-M = fit(ICA, X, k; do_whiten=false, tol=Inf)
-@test isa(M, ICA)
-@test indim(M) == m
-@test outdim(M) == k
-@test mean(M) == mv
-W = M.W
-@test transform(M, X) ≈ W' * (X .- mv)
-@test W'W ≈ Matrix(I, k, k)
-
-M = fit(ICA, X, k; do_whiten=true, tol=Inf)
-@test isa(M, ICA)
-@test indim(M) == m
-@test outdim(M) == k
-@test mean(M) == mv
-W = M.W
-@test W'C * W ≈ Matrix(I, k, k)
-
-@test_throws StatsBase.ConvergenceException fit(ICA, X, k; do_whiten=true, tol=1)
-
-# Use data of different type
-
-XX = convert(Matrix{Float32}, X)
-
-M = fit(ICA, XX, k; do_whiten=true, tol=Inf)
-@test eltype(mean(M)) == Float32
-@test eltype(M.W) == Float32
-
-M = fit(ICA, XX, k; do_whiten=false, tol=Inf)
-@test isa(M, ICA)
-@test eltype(mean(M)) == Float32
-@test eltype(M.W) == Float32
-W = M.W
-@test transform(M, X) ≈ W' * convert(Matrix{Float32}, (X .- mv))
-@test W'W ≈ Matrix{Float32}(I, k, k)
+@testset "ICA" begin
+
+    Random.seed!(15678)
+
+    function generatetestdata(n, k, m)
+        t = range(0.0, step=10.0, length=n)
+        s1 = sin.(t * 2)
+        s2 = s2 = 1.0 .- 2.0 * Bool[isodd(floor(Int, x / 3)) for x in t]
+        s3 = Float64[mod(x, 5.0) for x in t]
+
+        s1 += 0.1 * randn(n)
+        s2 += 0.1 * randn(n)
+        s3 += 0.1 * randn(n)
+
+        S = hcat(s1, s2, s3)'
+        @assert size(S) == (k, n)
+
+        A = randn(m, k)
+
+        X = A * S
+        mv = vec(mean(X, dims=2))
+        @assert size(X) == (m, n)
+        C = cov(X, dims=2)
+        return X, mv, C
+    end
+
+    @testset "Auxiliary" begin
+
+        f = icagfun(:tanh)
+        u, v = evaluate(f, 1.5)
+        @test u ≈ 0.905148253644866438242
+        @test v ≈ 0.180706638923648530597
+
+        f = icagfun(:tanh, Float32)
+        u, v = evaluate(f, 1.5f0)
+        @test u ≈ 0.90514827f0
+        @test v ≈ 0.18070662f0
+
+        f = icagfun(:tanh, 1.5)
+        u, v = evaluate(f, 1.2)
+        @test u ≈ 0.946806012846268289646
+        @test v ≈ 0.155337561057228069719
+
+        f = icagfun(:tanh, 1.5f0)
+        u, v = evaluate(f, 1.2f0)
+        @test u ≈ 0.94680610f0
+        @test v ≈ 0.15533754f0
+
+        f = icagfun(:gaus)
+        u, v = evaluate(f, 1.5)
+        @test u ≈ 0.486978701037524594696
+        @test v ≈ -0.405815584197937162246
+
+        f = icagfun(:gaus, Float32)
+        u, v = evaluate(f, 1.5f0)
+        @test u ≈ 0.4869787f0
+        @test v ≈ -0.40581557f0
+
+    end
+
+    ## data
+    @testset "Algorithm" begin
+
+        # sources
+        n = 1000
+        k = 3
+        m = 8
+        X, μ, C = generatetestdata(n, k, m)
+
+        # FastICA
+
+        M = fit(ICA, X, k; do_whiten=false, tol=Inf)
+        @test isa(M, ICA)
+        @test indim(M) == m
+        @test outdim(M) == k
+        @test mean(M) == μ
+        W = M.W
+        @test transform(M, X) ≈ W' * (X .- μ)
+        @test W'W ≈ Matrix(I, k, k)
+
+        M = fit(ICA, X, k; do_whiten=true, tol=Inf)
+        @test isa(M, ICA)
+        @test indim(M) == m
+        @test outdim(M) == k
+        @test mean(M) == μ
+        W = M.W
+        @test W'C * W ≈ Matrix(I, k, k)
+
+        @test_throws StatsBase.ConvergenceException fit(ICA, X, k; do_whiten=true, tol=0.01)
+
+        # Use data of different type
+
+        XX = convert(Matrix{Float32}, X)
+
+        M = fit(ICA, XX, k; do_whiten=true, tol=Inf)
+        @test eltype(mean(M)) == Float32
+        @test eltype(M.W) == Float32
+
+        M = fit(ICA, XX, k; do_whiten=false, tol=Inf)
+        @test isa(M, ICA)
+        @test eltype(mean(M)) == Float32
+        @test eltype(M.W) == Float32
+        W = M.W
+        @test transform(M, X) ≈ W' * convert(Matrix{Float32}, (X .- μ))
+        @test W'W ≈ Matrix{Float32}(I, k, k)
+
+        # input as view
+        M = fit(ICA, view(XX, :, 1:400), k; do_whiten=true, tol=Inf)
+        @test eltype(mean(M)) == Float32
+        @test eltype(M.W) == Float32
+    end
+
+end