-
Notifications
You must be signed in to change notification settings - Fork 7
/
nearcorr_aa.m
241 lines (209 loc) · 8.83 KB
/
nearcorr_aa.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
function [X,iter] = nearcorr_aa(A,pattern,mMax,itmax,ls_solve,delta,...
tol,droptol,AAstart)
%nearcorr_aa Nearest correlation matrix with Anderson Acceleration.
% [X,ITER] = nearcorr_aa(A,PATTERN,MMAX,ITMAX,LS_SOLVE,DELTA,...
% TOL,DROPTOL,BETA,AASTART)
% finds the nearest correlation matrix X to a symmetric matrix A by
% the alternating projections method with Anderson acceleration.
% ITER is the number of iterations taken.
% Prescribed elements can be kept fixed by specifying them in the
% matrix PATTERN: the (i,j) element is 1 if the corresponding element
% of A is to remain fixed, else zero.
% By default PATTERN is empty, meaning no elements are fixed.
% If PATTERN is non-empty then the unit diagonal must be explicitly
% forced if required.
% The solution can be forced to be positive definite with smallest
% eigenvalue at least DELTA, 0 < DELTA <= 1. Default: DELTA = 0.
% MMAX = history length parameter (non-negative integer).
% Default: 2. MMAX = 0 means no acceleration.
% ITMAX = maximum allowable number of iterations. Default: 100.
% LS_SOLVE = the least-squares solve method:
% 'u': QR factorization with updating (default),
% 'n': normal equations,
% 'b': MATLAB's backslash.
% TOL = convergence tolerance. Default: length(A)*eps.
% DROPTOL = tolerance for dropping stored residual vectors to improve
% conditioning. If DROPTOL > 0, drop residuals if the
% condition number exceeds DROPTOL; if DROPTOL <= 0,
% do not drop residuals. Default: 0.
% AASTART = acceleration delay factor: If AASTART > 0, start
% acceleration when iter = AAstart. Default: 1.
% Based on the code AndAcc.m by Homer Walker dated 10/14/2011,
% which is listed in
% Homer F. Walker, Anderson acceleration: Algorithms and implementations.
% Technical Report MS-6-15-50, Mathematical Sciences Department,
% Worcester Polytechnic Institute, Worcester, MA 01609, USA, June 2011.
% This code is specialized for the alternating projections method for the
% nearest correlation matrix, adds two extra options for solving
% the least squares problem, and contains many other changes.
% Nick Higham and Natasa Strabic, 2015.
if ~isequal(A,A'), error('The input matrix must by symmetric.'), end
n = length(A);
if nargin < 2, pattern = []; end
if nargin < 3 || isempty(mMax), mMax = 2; end
if nargin < 4 || isempty(itmax), itmax = 100; end
if nargin < 5 || isempty(ls_solve), ls_solve = 'u'; end
if nargin < 6 || isempty(delta), delta = 0; end
if nargin < 7 || isempty(tol), tol = n*eps; end
if nargin < 8 || isempty(droptol), droptol = 0; end
if nargin < 9, AAstart = 1; end
% Initialize the storage arrays.
DG = []; % Storage of g-value differences.
DF = []; % Storage of f-value differences, needed for 'n' and 'b' options.
% Initialize the number of stored residuals.
mAA = 0;
% Initialization for the nearest corr. matrix problem.
Yin = A;
Sin = zeros(n);
% Top of the iteration loop.
for iter = 1:itmax
% Apply g
[Xout,Yout,Sout] = ap_step(A,Yin,Sin,pattern,delta);
% We have to vectorize the output.
x = [Yin(:); Sin(:)];
gval = [Yout(:); Sout(:)];
fval = gval - x;
% Convergence test.
rel_diffXY = norm(Yout - Xout,'fro')/norm(Yout,'fro');
if rel_diffXY < tol, break; end
if mMax == 0 || iter < AAstart,
% Without acceleration, update x <- g(x) to obtain the next
% approximate solution.
x = gval;
else
% Apply Anderson acceleration.
% Update the df vector and the DG array.
% If the least-squares solves are by normal equations or backslash, we
% need DF, the matrix of function value differences built in the same
% way as DG.
if iter > AAstart
df = fval-f_old;
if mAA < mMax
DG = [DG gval-g_old];
if ls_solve ~= 'u'
DF = [DF df];
end
else
DG(:,1)=[];
DG(:,mAA) = gval-g_old;
if ls_solve ~= 'u'
DF(:,1) = [];
DF(:,mAA) = df;
end
end
mAA = mAA + 1;
end
f_old = fval;
g_old = gval;
if mAA == 0
% Update x <- g(x) to obtain the next approximate
% solution. This is the x1 = g(x0) line in pseudocode.
x = gval;
else
% Solve least-squares problem and update solution.
switch ls_solve
case 'u' % Use the QR factors.
if mAA == 1
% Form the initial QR decomposition.
R(1,1) = norm(df);
Q = R(1,1)\df;
else
% Update the QR decomposition.
if mAA > mMax
% If the column dimension of Q is mMax, delete
% the first column and update the decomposition.
[Q,R] = qrdelete(Q,R,1);
mAA = mAA - 1;
% The following treats the qrdelete quirk described below.
if size(R,1) ~= size(R,2),
Q = Q(:,1:mAA-1); R = R(1:mAA-1,:);
end
% Explanation: If Q is not square, then
% qrdelete(Q,R,1) reduces the column
% dimension of Q by 1 and the column and row
% dimensions of R by 1. But if Q *is* square, then
% the column dimension of Q is not reduced and only
% the column dimension of R is reduced by one. This
% is to allow for MATLAB’s default "thick" QR
% decomposition, which always produces a square Q.
end
% Now update the QR decomposition to incorporate
% the new column.
for j = 1:mAA - 1
R(j,mAA) = Q(:,j)'*df;
df = df - R(j,mAA)*Q(:,j);
end
R(mAA,mAA) = norm(df);
Q = [Q,R(mAA,mAA)\df];
end
if droptol > 0
% Drop residuals to improve conditioning if necessary.
condDF = cond(R);
while condDF > droptol && mAA > 1
fprintf('cond(D) = %e, reducing mAA to %d \n', ...
condDF, mAA-1);
[Q,R] = qrdelete(Q,R,1);
DG = DG(:,2:mAA);
mAA = mAA - 1;
% Treat the qrdelete quirk described above.
if size(R,1) ~= size(R,2),
Q = Q(:,1:mAA); R = R(1:mAA,:);
end
condDF = cond(R);
end
end % if droptolerance
% Solve the least-squares problem.
gamma = R\(Q'*fval);
case 'n' % use normal equations
if mAA > mMax, mAA = mAA - 1; end
c = DF'*fval;
[~,R] = qr(DF,0); % economy size QR
y = R'\c;
gamma = R\y;
case 'b' % use backslash
if mAA > mMax, mAA = mAA - 1; end
gamma = DF\fval;
end % switch
% Update the approximate solution.
x = gval - DG*gamma;
end % end of least-squares solve, we have the next iterate x
end % end of the acceleration part
% Update the parameters: undo reshape of x to get new Yin, Sin.
Yin(:) = x(1:n^2);
Sin(:) = x(n*n+1:2*n^2);
end
% Bottom of the iteration loop.
if rel_diffXY > tol && iter == itmax,
error(['Stopped after ' num2str(itmax) ' its. Try increasing ITMAX.'])
end
X = Yin;
end
% --------------------------------------------
% Subfunctions
function [Xout,Yout,Sout] = ap_step(A,Yin,Sin,pattern,delta)
% Fixed point iteration for the alternating projections step
% for the nearest correlation matrix.
R = Yin - Sin;
Xout = proj_spd(R,delta);
Sout = Xout - R;
Yout = proj_pattern(A,Xout,pattern);
end
function X = proj_pattern(A,X,pattern)
% Return the nearest matrix to X with fixed elements from A specified in
% the positions PATTERN
n = length(A);
if isempty(pattern)
% Standard NCM, with no fixed elements. Set unit diagonal.
X(1:n+1:n^2) = 1;
else
F = find(pattern); % Locations of fixed elements in A.
X(F) = A(F);
end
end
function X = proj_spd(A,delta)
% Return the nearest positive semidefinite matrix to A with smallest
% eigenvalue at least delta.
[V,D] = eig(A);
X = V*diag(max(diag(D),delta))*V';
X = (X+X')/2; % Ensure symmetry.
end