-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathRPADMMQP.m
64 lines (64 loc) · 1.43 KB
/
RPADMMQP.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
%
function [x,y]=RPADMMQP(Q,A,b,c,nb)
tic
toler=1.e-4;
[m,n]=size(A);
beta=1;
x=ones(n,1);
s=x;
y=zeros(m,1);
z=zeros(n,1);
nv=floor(n/nb);
ATA=zeros(nv,n);
% store block-wise inverse
for ii=1:nb,
il=(ii-1)*nv+1;
iu=il+nv-1;
AB=A(:,il:iu);
ATA(:,il:iu)=inv(Q(il:iu,il:iu)+beta*(AB'*AB)+beta*eye(nv));
end;
clear AB;
% greate two vectos: the gradient of objective and constraint LHS
Qxc=Q*x+c;
Axb=A*x-b;
%
iter=0;
for k=1:300,
or=randperm(nb)
% Update x following a random order
for ii=1:nb,
il=(or(ii)-1)*nv+1
iu=il+nv-1
Qxc=Qxc-Q(:,il:iu)*x(il:iu);
Axb=Axb-A(:,il:iu)*x(il:iu);
cc=Qxc(il:iu)-z(il:iu)-beta*s(il:iu)-A(:,il:iu)'*(y-beta*Axb);
x(il:iu)=-ATA(:,il:iu)*cc;
Qxc=Qxc+Q(:,il:iu)*x(il:iu);
Axb=Axb+A(:,il:iu)*x(il:iu);
end
% Update s
cc=x-(1/beta)*z;
s=max(0,cc);
% Update multipliers
y=y-beta*Axb;
z=z-beta*(x-s);
%
iter=iter+1;
% Check stopping criteria
if (norm(Axb)+norm(x-s))/(1+norm(x))<toler, break, end;
end;
toc
x=s;
iter=iter
% Randomly Permuted Multi-Block ADMM for solving convex quadratic program
%
% minimize 0.5x'Qx + c'x
% subject to Ax = b, (y dimension m)
% x - s = 0, (z dimesnion n)
% s>=0
%
% Input: A, Q, b, c, nb(number of blocks such that number of variables
% divided by nb is an integer)
%
% Output:primal x, dual y
%