1function [IR, CC, D] = adim(U, UC, IR, CC, arg5);
2% ADIM adaptive information matrix. Estimates the inverse
3%   correlation matrix in an adaptive way.
4%
5% [IR, CC] = adim(U, UC [, IR0 [, CC0]]);
6%   U 	input data
7%   UC 	update coefficient 0 < UC << 1
8%   IR0	initial information matrix
9%   CC0 initial correlation matrix
10%   IR	information matrix (inverse correlation matrix)
11%   CC  correlation matrix
12%
13%  The algorithm uses the Matrix Inversion Lemma, also known as
14%     "Woodbury's identity", to obtain a recursive algorithm.
15%     IR*CC - UC*I should be approx. zero.
16%
17% Reference(s):
18% [1] S. Haykin. Adaptive Filter Theory, Prentice Hall, Upper Saddle River, NJ, USA
19%     pp. 565-567, Equ. (13.16), 1996.
20
21
22%       $Id$
23%       Copyright (C) 1998-2003 by Alois Schloegl <alois.schloegl@gmail.com>
24%
25%    This program is free software: you can redistribute it and/or modify
26%    it under the terms of the GNU General Public License as published by
27%    the Free Software Foundation, either version 3 of the License, or
28%    (at your option) any later version.
29%
30%    This program is distributed in the hope that it will be useful,
31%    but WITHOUT ANY WARRANTY; without even the implied warranty of
32%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
33%    GNU General Public License for more details.
34%
35%    You should have received a copy of the GNU General Public License
36%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
37
38
39
40[ur,p] = size(U);
41
42Mode_E = 1;
43%if nargin < 4,
44%        m  = zeros(size(U)+[0,Mode_E]);
45%end;
46
47if nargin<2,
48        fprintf(2,'Error ADIM: missing update coefficient\n');
49        return;
50else
51        if ~((UC > 0) & (UC <1)),
52                fprintf(2,'Error ADIM: update coefficient not within range [0,1]\n');
53                return;
54        end;
55        if UC > 1/p,
56                fprintf(2,'Warning ADIM: update coefficient should be smaller than 1/number_of_dimensions\n');
57        end;
58end;
59
60if nargin<3,
61        IR = [];
62end;
63if nargin<4,
64        CC = [];
65end;
66if nargin<5,
67        arg5 = 6;
68end;
69if isempty(IR),
70        IR = eye(p+Mode_E);
71end;
72if isempty(CC),
73        CC = eye(p+Mode_E);
74end;
75
76D = zeros(ur,(p+Mode_E)^2);
77%D = zeros(ur,1);
78W  = eye(p+Mode_E)*UC/(p+Mode_E);
79W2 = eye(p+Mode_E)*UC*UC/(p+Mode_E);
80
81for k = 1:ur,
82        if ~Mode_E,
83                % w/o mean term
84                d  = U(k,:);
85        else
86                % w/ mean term
87                d  = [1,U(k,:)];
88        end;
89
90        if ~any(isnan(d)),
91                CC = (1-UC)*CC + UC*(d'*d);
92
93                %K = (1+UC)*IR*d'/(1+(1+UC)*d*IR*d');
94                v  = IR*d';
95                %K = v/(1-UC+d*v);
96
97                if arg5==0;  		% original version according to [1], can become unstable
98                        IR = (1+UC)*IR - (1+UC)/(1-UC+d*v)*v*v';
99
100                elseif arg5==1;  	% this provides IR*CC == I, this seems to be stable
101                        IR = IR - (1+UC)/(1-UC+d*v)*v*v' + sum(diag(IR))*W;
102
103                elseif arg5==6;  	% DEFAULT:
104                        IR = ((1+UC)/2)*(IR+IR') - ((1+UC)/(1-UC+d*v))*v*v';
105
106                        % more variants just for testing, do not use them.
107                elseif arg5==2;
108                        IR = IR - (1+UC)/(1-UC+d*v)*v*v' + sum(diag(IR))*W2;
109
110                elseif arg5==3;
111                        IR = IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E);
112
113                elseif arg5==4;
114                        IR = (1+UC)*IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E);
115
116                elseif arg5==5;
117                        IR = IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E);
118
119                end;
120        end;
121
122        %D(k) = det(IR);
123        D(k,:) = IR(:)';
124end;
125
126IR = IR / UC;
127if Mode_E,
128        IR(1,1) = IR(1,1) - 1;
129end;
130
131
132