Skip to content

Commit

Permalink
ICA: Fixed convergence criteria (closes #45)
Browse files Browse the repository at this point in the history 
- Added testsets.
- Removed `Printf` dependency.
  • Loading branch information
@wildart
wildart authored Oct 2, 2018
1 parent 4b6b8d2 commit 2c61429
Show file tree
Hide file tree
Showing 2 changed files with 123 additions and 127 deletions.
40 changes: 10 additions & 30 deletions src/ica.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand Down
210 changes: 113 additions & 97 deletions test/ica.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit 2c61429

Please sign in to comment.