-
Notifications
You must be signed in to change notification settings - Fork 3
/
icaML.m
783 lines (355 loc) · 15.1 KB
/
icaML.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
function [S,A,U,ll,info]=icaML(X,K,par,debug_draw)
% icaML : ICA by ML (Infomax) with square mixing matrix and no noise.
%
% function [S,A,U,ll,info]=icaML(X,[K],[par]) Independent component analysis (ICA) using
% maximum likelihood, square mixing matrix and
% no noise [1] (Infomax). Source prior is assumed
% to be p(s)=1/pi*exp(-ln(cosh(s))). For optimization
% the BFGS algorithm is used [2]. See code for
% references.
%
% X : Mixed signals
% K : Number of source components.
% For K=0 (default) number of sources are equal to number of
% observations.
% For K < number of observations, SVD is used to reduce the
% dimension.
% par: Vector with 4 elements:
% (1) : Expected length of initial step
% Stopping criteria:
% (2) : Gradient ||g||_inf <= par(2)
% (3) : Parameter changes ||dW||_2 <= par(3)*(par(3) + ||W||_2)
% (4) : Maximum number of iterations
% Any illegal element in opts is replaced by its
% default value, [1 1e-4*||g(x0)||_inf 1e-8 100]
% debug_draw : Draw debug information
%
% S : Estimated source signals with variance
% scaled to one.
% A : Estimated mixing matrix
% U : Principal directions of preprocessing PCA.
% If K (the number of sources) is equal to the number
% of observations then no PCA is performed and U=eye(K).
% ll : Log likelihood for estimated sources
% info : Performance information, vector with 6 elements:
% (1:3) : final values of [ll ||g||_inf ||dx||_2]
% (4:5) : no. of iteration steps and evaluations of (ll,g)
% (6) : 1 means stopped by small gradient
% 2 means stopped by small x-step
% 3 means stopped by max number of iterations.
% 4 means stopped by zero step.
%
%
% Eks. Separate with number of sources equal number of
% observations.
%
% [S,A] = icaML(X);
%
% Eks. Separate with number of sources equal to k, using SVD
% as pre-processing.
%
% [S,A,U] = icaML(X,k);
%
% - version 1.5 (Revised 9/9-2003)
% - IMM, Technical University of Denmark
% - version 1.4
% - by Thomas Kolenda 2002 - IMM, Technical University of Denmark
% Revised: 9/9-2003, Mads, [email protected]
% * Fixed help message to inform the user about U.
% * Removed the automatic use of PCA in the quadratic case.
% Bibtex references:
% [1]
% @article{Bell95,
% author = "A. Bell and T.J. Sejnowski",
% title = "An Information-Maximization Approach to Blind Separation and Blind Deconvolution",
% journal = "Neural Computation",
% year = "1995",
% volume = "7",
% pages = "1129-1159",
% }
%
% [2]
% @techreport{Nielsen01:unopt,
% author = "H.B. Nielsen",
% title = "UCMINF - an Algorithm for Unconstrained, Nonlinear Optimization ",
% institution = "IMM, Technical University of Denmark",
% number = "IMM-TEC-0019",
% year = "2001",
% url = "http://www.imm.dtu.dk/pubdb/views/edoc_download.php/642/ps/imm642.ps",
% }
% algorithm parameter settings
MaxNrIte = 1000;
try
debug = debug_draw;
catch
debug = 0;
end
if debug==1,fprintf('\n** Start icaML ***************************************\n');tic;end
% Scale X to avoid numerical problems
Xorg = X;
scaleX = max([abs(max(X(:))),abs(min(X(:)))]);
X = X./scaleX;
% set number of source parameters
if nargin<2, K=size(X,1); end
if ((K>0) & (K<size(X,1))),
[U,X] = callSVD(X,K,debug);
else
U=eye(size(X,1));
if debug==1,disp('Don''t use SVD');end
end
% initialize optimize parameters
try
ucminf_opt(1) = par(1);
ucminf_opt(2) = par(2);
ucminf_opt(3) = par(3);
ucminf_opt(4) = par(4);
catch
if debug==1,disp('Use default parameters to optimize');end
ucminf_opt = [1 1e-4 1e-8 MaxNrIte];
end
% initialize variables
[M,N] = size(X);
W = eye(M);%W= rand(M);
if debug==1,fprintf('Number of samples %d - Number of sources %d\n',N,M);end
% optimize
if debug==1,fprintf('Optimize ICA ... ');end
par.X=X; par.M=M; par.N=N;
[W,info] = ucminf( 'ica_MLf' , par , W(:) , ucminf_opt );
W = reshape(W,M,M);
% estimates
A=pinv(W);
S=W*X;
if debug==1,disp('done optimize ICA!');end
% sort components according to energy
Avar=diag(A'*A)/M;
Svar=diag(S*S')/N;
vS=var(S');
sig=Avar.*Svar;
[a,indx]=sort(sig);
S=S(indx(M:-1:1),:);
A=A(:,indx(M:-1:1));
% scale back
A=A.*scaleX;
% log likelihood
if nargout>3
ll= N*log(abs(det(inv(A)))) - sum(sum(log( cosh(S) ))) - N*M*log(pi);
end
if debug==1,fprintf('** End of icaML ****** time %0.2f sec******************\n\n',toc/60);end
function [f,dW]=ica_MLf(W,par)
% returns the negative log likelihood and its gradient w.r.t. W
X=par.X; M=par.M; N=par.N;
W=reshape(W,M,M);
S=W*X;
% negative log likelihood function
f=-( N*log(abs(det(W))) - sum(sum(log( cosh(S) ))) - N*M*log(pi) );
if nargout>1
% gradient w.r.t. W
dW=-(N*inv(W') - tanh(S)*X');
dW=dW(:);
end
function [U,DV] = callSVD(X,K,draw)
% Reduce dimension with SVD
[M,N] = size(X);
if N>M % Transpose if matrix is flat, to speed up the svd and later transpose back again
if draw==1,disp('Do Transpose SVD');end
[V,D,U] = svd(X',0);
else;
if draw==1,disp('Do SVD');end
[U,D,V] = svd(X,0);
end;
DV = D(1:K,1:K)*V(:,1:K)';
function [X, info, perf, D] = ucminf(fun,par, x0, opts, D0)
%UCMINF BFGS method for unconstrained nonlinear optimization:
% Find xm = argmin{f(x)} , where x is an n-vector and the scalar
% function F with gradient g (with elements g(i) = DF/Dx_i )
% must be given by a MATLAB function with with declaration
% function [F, g] = fun(x, par)
% par holds parameters of the function. It may be dummy.
%
% Call: [X, info {, perf {, D}}] = ucminf(fun,par, x0, opts {,D0})
%
% Input parameters
% fun : String with the name of the function.
% par : Parameters of the function. May be empty.
% x0 : Starting guess for x .
% opts : Vector with 4 elements:
% opts(1) : Expected length of initial step
% opts(2:4) used in stopping criteria:
% ||g||_inf <= opts(2) or
% ||dx||_2 <= opts(3)*(opts(3) + ||x||_2) or
% no. of function evaluations exceeds opts(4) .
% Any illegal element in opts is replaced by its
% default value, [1 1e-4*||g(x0)||_inf 1e-8 100]
% D0 : (optional) If present, then approximate inverse Hessian
% at x0 . Otherwise, D0 := I
% Output parameters
% X : If perf is present, then array, holding the iterates
% columnwise. Otherwise, computed solution vector.
% info : Performance information, vector with 6 elements:
% info(1:3) = final values of [f(x) ||g||_inf ||dx||_2]
% info(4:5) = no. of iteration steps and evaluations of (F,g)
% info(6) = 1 : Stopped by small gradient
% 2 : Stopped by small x-step
% 3 : Stopped by opts(4) .
% 4 : Stopped by zero step.
% perf : (optional). If present, then array, holding
% perf(1:2,:) = values of f(x) and ||g||_inf
% perf(3:5,:) = Line search info: values of
% alpha, phi'(alpha), no. fct. evals.
% perf(6,:) = trust region radius.
% D : (optional). If present, then array holding the
% approximate inverse Hessian at X(:,end) .
% Hans Bruun Nielsen, IMM, DTU. 00.12.18 / 02.01.22
% Check call
[x n F g] = check(fun,par,x0,opts);
if nargin > 4, D = checkD(n,D0); fst = 0;
else, D = eye(n); fst = 1; end
% Finish initialization
k = 1; kmax = opts(4); neval = 1; ng = norm(g,inf);
Delta = opts(1);
Trace = nargout > 2;
if Trace
X = x(:)*ones(1,kmax+1);
perf = [F; ng; zeros(3,1); Delta]*ones(1,kmax+1);
end
found = ng <= opts(2);
h = zeros(size(x)); nh = 0;
ngs = ng * ones(1,3);
while ~found
% Previous values
xp = x; gp = g; Fp = F; nx = norm(x);
ngs = [ngs(2:3) ng];
h = D*(-g(:)); nh = norm(h); red = 0;
if nh <= opts(3)*(opts(3) + nx), found = 2;
else
if fst | nh > Delta % Scale to ||h|| = Delta
h = (Delta / nh) * h; nh = Delta;
fst = 0; red = 1;
end
k = k+1;
% Line search
[al F g dval slrat] = softline(fun,par,x,F,g, h);
if al < 1 % Reduce Delta
Delta = .35 * Delta;
elseif red & (slrat > .7) % Increase Delta
Delta = 3*Delta;
end
% Update x, neval and ||g||
x = x + al*h; neval = neval + dval; ng = norm(g,inf);
if Trace
X(:,k) = x(:);
perf(:,k) = [F; ng; al; dot(h,g); dval; Delta]; end
h = x - xp; nh = norm(h);
if nh == 0,
found = 4;
else
y = g - gp; yh = dot(y,h);
if yh > sqrt(eps) * nh * norm(y)
% Update D
v = D*y(:); yv = dot(y,v);
a = (1 + yv/yh)/yh; w = (a/2)*h(:) - v/yh;
D = D + w*h' + h*w';
end % update D
% Check stopping criteria
thrx = opts(3)*(opts(3) + norm(x));
if ng <= opts(2), found = 1;
elseif nh <= thrx, found = 2;
elseif neval >= kmax, found = 3;
% elseif neval > 20 & ng > 1.1*max(ngs), found = 5;
else, Delta = max(Delta, 2*thrx); end
end
end % Nonzero h
end % iteration
% Set return values
if Trace
X = X(:,1:k); perf = perf(:,1:k);
else, X = x; end
info = [F ng nh k-1 neval found];
% ========== auxiliary functions =================================
function [x,n, F,g, opts] = check(fun,par,x0,opts0)
% Check function call
x = x0(:); sx = size(x); n = max(sx);
if (min(sx) > 1)
error('x0 should be a vector'), end
[F g] = feval(fun,x,par);
sf = size(F); sg = size(g);
if any(sf-1) | ~isreal(F)
error('F should be a real valued scalar'), end
if (min(sg) ~= 1) | (max(sg) ~= n)
error('g should be a vector of the same length as x'), end
so = size(opts0);
if (min(so) ~= 1) | (max(so) < 4) | any(~isreal(opts0(1:4)))
error('opts should be a real valued vector of length 4'), end
opts = opts0(1:4); opts = opts(:).';
i = find(opts <= 0);
if length(i) % Set default values
d = [1 1e-4*norm(g, inf) 1e-8 100];
opts(i) = d(i);
end
% ---------- end of check ---------------------------------------
function D = checkD(n,D0)
% Check given inverse Hessian
D = D0; sD = size(D);
if any(sD - n)
error(sprintf('D should be a square matrix of size %g',n)), end
% Check symmetry
dD = D - D'; ndD = norm(dD(:),inf);
if ndD > 10*eps*norm(D(:),inf)
error('The given D0 is not symmetric'), end
if ndD, D = (D + D')/2; end % Symmetrize
[R p] = chol(D);
if p
error('The given D0 is not positive definite'), end
function [alpha,fn,gn,neval,slrat] = ... softline(fun,fpar, x,f,g, h)
% Soft line search: Find alpha = argmin_a{f(x+a*h)}
% Default return values
alpha = 0; fn = f; gn = g; neval = 0; slrat = 1;
n = length(x);
% Initial values
dfi0 = dot(h,gn); if dfi0 >= 0, return, end
fi0 = f; slope0 = .05*dfi0; slopethr = .995*dfi0;
dfia = dfi0; stop = 0; ok = 0; neval = 0; b = 1;
while ~stop
[fib g] = feval(fun,x+b*h,fpar); neval = neval + 1;
dfib = dot(g,h);
if b == 1, slrat = dfib/dfi0; end
if fib <= fi0 + slope0*b % New lower bound
if dfib > abs(slopethr), stop = 1;
else
alpha = b; fn = fib; gn = g; dfia = dfib;
ok = 1; slrat = dfib/dfi0;
if (neval < 5) & (b < 2) & (dfib < slopethr)
% Augment right hand end
b = 2*b;
else, stop = 1; end
end
else, stop = 1; end
end
stop = ok; xfd = [alpha fn dfia; b fib dfib; b fib dfib];
while ~stop
c = interpolate(xfd,n);
[fic g] = feval(fun, x+c*h, fpar); neval = neval+1;
xfd(3,:) = [c fic dot(g,h)];
if fic < fi0 + slope0*c % New lower bound
xfd(1,:) = xfd(3,:); ok = 1;
alpha = c; fn = fic; gn = g; slrat = xfd(3,3)/dfi0;
else, xfd(2,:) = xfd(3,:); ok = 0; end
% Check stopping criteria
ok = ok & abs(xfd(3,3)) <= abs(slopethr);
stop = ok | neval >= 5 | diff(xfd(1:2,1)) <= 0;
end % while
%------------ end of softline ------------------------------
function alpha = interpolate(xfd,n);
% Minimizer of parabola given by
% xfd(1:2,1:3) = [a fi(a) fi'(a); b fi(b) dummy]
a = xfd(1,1); b = xfd(2,1); d = b - a; dfia = xfd(1,3);
C = diff(xfd(1:2,2)) - d*dfia;
if C >= 5*n*eps*b % Minimizer exists
A = a - .5*dfia*(d^2/C);
d = 0.1*d; alpha = min(max(a+d, A), b-d);
else
alpha = (a+b)/2;
end
%------------ end of interpolate --------------------------