-
Notifications
You must be signed in to change notification settings - Fork 6
/
blockCoordinateDescent.m
51 lines (49 loc) · 1.83 KB
/
blockCoordinateDescent.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
function [x, multiplier] = blockCoordinateDescent(c, mm, a, lambda, reserve)
M = size(mm, 1);
if norm(c, 2) <= lambda / 2 || a < 0.9 * reserve
x = zeros(M, 1);
multiplier = 0;
else
theta = 1 / sqrt(a);
multiplier = 0;
target = bisectionTarget(0, theta, lambda, mm, c);
if target > 1
miuLow = 0;
miuHigh = norm(c) * theta * 1.1;
multiplier = (miuLow + miuHigh) / 2;
theta = 1 / sqrt(a);
target = bisectionTarget(multiplier, theta, lambda, mm, c);
while abs(target - 1) > 1e-6
if target > 1
miuLow = multiplier;
elseif target < 1
miuHigh = multiplier;
end
multiplier = (miuLow + miuHigh) / 2;
target = bisectionTarget(multiplier, theta, lambda, mm, c);
end
elseif target < 1
tLow = 0;
tHigh = (spectralRadius(mm) + norm(c, 2) * theta * 1.1) / (norm(c, 2) - lambda / 2);
multiplier = 0;
theta = (tLow + tHigh) / 2;
target = bisectionTarget(multiplier, theta, lambda, mm, c);
while abs(target - 1) > 1e-6
if target > 1
tHigh = theta;
elseif target < 1
tLow = theta;
end
theta = (tLow + tHigh) / 2;
target = bisectionTarget(multiplier, theta, lambda, mm, c);
end
end
x = (mm + (lambda * theta / 2 + multiplier) * eye(size(mm))) \ c;
end
return
function rho = spectralRadius(A)
rho = max(abs(eig(A)));
return
function r = bisectionTarget(multiplier, theta, lambda, mm, c)
r = theta * norm((mm + (lambda * theta / 2 + multiplier) * eye(size(mm))) \ c);
return