-
Notifications
You must be signed in to change notification settings - Fork 0
/
rmlr.m
152 lines (126 loc) · 4.36 KB
/
rmlr.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
% Name : Robust Multinomial Logistic Regression
% Author: Jakramate Bootkrajang
% Last update: 16 April 2013
% Input :
% w = given initial model's parameters (DIM,CLS)
% g = initial label flipping probability matrix
% x = Design matrix where a row represents a sample
% y = Target value
% options.maxIter = maximum optimising iteration
% options.estGam = estimating the gamma using multiplicative update
% options.regFunc = type of regularisation: 'noreg', 'L1' or 'L2'
% options.sn = Small number for approximating L1 objective.
% options.verbose = displaying negative log-likelihood
% Output:
% w = fitted model
% g = estimated gamma matrix
% llh = negative log-likelihood
% Note :
% 1. The model uses {1,2,...K} class representation internally.
% |-----|-----------|
% | | 1 y^ n|
% ------|-----------|
% |y 1 | g00 g0n | where n = k-1 and k = number of class
% | n | gn0 gnn |
% |-----|-----------|
%=====================================================================================
function [w, g, llh] = rmlr(w, g, x, y, options)
% check if bias term is added i.e. first
% column must be unique and equal to 1
if (sum(x(:,1)) ~= size(x,1))
disp('Bias terms might have not been added');
end
if ~isfield(options,'sn')
options.sn = 1e-8;
end
if ~isfield(options,'maxIter')
options.maxIter = 50;
end
CLS = length(unique(y));
DIM = size(x,2);
% create label mask
mask = repmat(y, 1, CLS);
for k = 1:CLS
mask(:,k) = logical(y==k);
end
mask = logical(mask);
alpha = zeros(size(g)); % Bayesian approach
% ========================= BEGIN ESTIMATING PARAMETERS ===========================
for l=1:options.maxIter
switch options.regFunc
case 'lasso'
lambda = 1 ./(sqrt(w.^2 + options.sn));
lambda(1,:) = 0;
case 'l2'
lambda = repmat((length(w)./2 +1)./(sum(w(2:end,:).^2)./2 + 2), length(w), 1);
lambda(1,:) = 0;
case 'noreg'
lambda = zeros(size(w));
otherwise
disp('Invalid regularisation function. Now using L1');
options.regFunc = 'L1';
end
if l==1
lambda = lambda * 0;
end
% minimise objective function
w0 = cat(2, w(:));
lambda0 = cat(2, lambda(:));
opts.maxIter = 10;
opts.Display = false;
[w0, fv, v1, v2] = minFunc(@gradws, w0, opts, g, x, y, mask, lambda0, options);
w = reshape(w0, DIM, CLS);
% Bayesian approach
if options.breg
alpha = 2./(g + 1e-1);
end
%alpha
if options.estG
for gl=1:1
% update the gamma
sk = posteriorYHat(x, w, g);
lj = softmax(x * w);
num = zeros(CLS,CLS);
for k = 1:CLS
skc = repmat(sk(:,k), 1, CLS);
tmp = lj ./ skc;
% select y_hat_i = k with mask
tmp2 = tmp(mask(:,k),:);
num(:,k) = sum(tmp2,1) ;
if options.breg
num(:,k) = num(:,k) + alpha(:,k) ;
end
end
% calculate denominator
if options.breg
dnt = sum((num .* g), 2) ;
else
dnt = sum((num .* g), 2);
end
denom = repmat(dnt, 1, CLS);
% now last step for gamma
g = g .* (num ./ denom);
%g(g<0) = 0;
end
end
% save the progress
llh(l) = fv(end);
if options.verbose
disp(llh(l));
end
end % end optimisation loop
% ============================== END ESTIMATING PARAMETER ============================
% class identity assignment, solved with Hungarian algorithm
task = repmat(max(max(g)),CLS,CLS) - g;
[asm n] = munkres(task);
s = 1;
for loop = 1:CLS
[nu i] = max(asm(s,:));
if (asm(s,i) > asm(i,i))
asm = rowswap(asm, s, i);
g = rowswap(g, s, i);
w = rowswap(w',s, i)';
else
s = s + 1;
end
end