1 | % * This code was used in the following articles:
|
---|
2 | % * [1] Learning 3-D Scene Structure from a Single Still Image,
|
---|
3 | % * Ashutosh Saxena, Min Sun, Andrew Y. Ng,
|
---|
4 | % * In ICCV workshop on 3D Representation for Recognition (3dRR-07), 2007.
|
---|
5 | % * (best paper)
|
---|
6 | % * [2] 3-D Reconstruction from Sparse Views using Monocular Vision,
|
---|
7 | % * Ashutosh Saxena, Min Sun, Andrew Y. Ng,
|
---|
8 | % * In ICCV workshop on Virtual Representations and Modeling
|
---|
9 | % * of Large-scale environments (VRML), 2007.
|
---|
10 | % * [3] 3-D Depth Reconstruction from a Single Still Image,
|
---|
11 | % * Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng.
|
---|
12 | % * International Journal of Computer Vision (IJCV), Aug 2007.
|
---|
13 | % * [6] Learning Depth from Single Monocular Images,
|
---|
14 | % * Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng.
|
---|
15 | % * In Neural Information Processing Systems (NIPS) 18, 2005.
|
---|
16 | % *
|
---|
17 | % * These articles are available at:
|
---|
18 | % * http://make3d.stanford.edu/publications
|
---|
19 | % *
|
---|
20 | % * We request that you cite the papers [1], [3] and [6] in any of
|
---|
21 | % * your reports that uses this code.
|
---|
22 | % * Further, if you use the code in image3dstiching/ (multiple image version),
|
---|
23 | % * then please cite [2].
|
---|
24 | % *
|
---|
25 | % * If you use the code in third_party/, then PLEASE CITE and follow the
|
---|
26 | % * LICENSE OF THE CORRESPONDING THIRD PARTY CODE.
|
---|
27 | % *
|
---|
28 | % * Finally, this code is for non-commercial use only. For further
|
---|
29 | % * information and to obtain a copy of the license, see
|
---|
30 | % *
|
---|
31 | % * http://make3d.stanford.edu/publications/code
|
---|
32 | % *
|
---|
33 | % * Also, the software distributed under the License is distributed on an
|
---|
34 | % * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
|
---|
35 | % * express or implied. See the License for the specific language governing
|
---|
36 | % * permissions and limitations under the License.
|
---|
37 | % *
|
---|
38 | % */
|
---|
39 | function [w, alfa, status, history, T_nt_hist] = ... |
---|
40 | SigmoidLogBarrierSolver( Para, t_0, alpha_0, w_0, method, ptol, pmaxi, VERBOSE) |
---|
41 | % All arguments are optional. You just need A, b (and S, q for inequality) (see INPUT below) |
---|
42 | % |
---|
43 | % l1 penalty approximate Problem with inequality constrain Solver |
---|
44 | % |
---|
45 | % exact problems form: |
---|
46 | % |
---|
47 | % minimize norm( Aw - b, 1) st, Sw+q<=0 |
---|
48 | % |
---|
49 | % sigmoid approximate form: |
---|
50 | % |
---|
51 | % minimize sum_i (1/alfa)*( log( 1+exp(-alfa*( A(i,:)*w - b(i)))) + ... |
---|
52 | % log( 1+exp(alfa*( A(i,:)*w - b(i)))) ) |
---|
53 | % st, Sw+q <= 0 (log barrier approximation) |
---|
54 | % |
---|
55 | % where variable is w and problem data are A, b, S, and q. |
---|
56 | % |
---|
57 | % INPUT |
---|
58 | % |
---|
59 | % A : mxn matrix; |
---|
60 | % b : m vector; |
---|
61 | % S : lxn matrix; |
---|
62 | % q : l vector; |
---|
63 | % |
---|
64 | % method : string; search direction method type |
---|
65 | % 'cg' : conjugate gradients method, 'pcg' |
---|
66 | % 'pcg' : preconditioned conjugate gradients method |
---|
67 | % 'exact': exact method (default value) |
---|
68 | % ptol : scalar; pcg relative tolerance. if empty, use adaptive rule. |
---|
69 | % pmaxi : scalar: pcg maximum iteration. if empty, use default value (500).a |
---|
70 | % mu : approximate_closeness factor; bigger -> more accurate |
---|
71 | % VERBOSE : enable disp |
---|
72 | % |
---|
73 | % OUTPUT |
---|
74 | % |
---|
75 | % w : n vector; classifier |
---|
76 | % status : scalar; +1: success, -1: maxiter exceeded |
---|
77 | % history : |
---|
78 | % row 1) phi (objective value) |
---|
79 | % row 2) norm(gradient of phi) |
---|
80 | % row 3) cumulative cg iterations |
---|
81 | % |
---|
82 | % USAGE EXAMPLE |
---|
83 | % |
---|
84 | % [w,status] = (A, b, S, q, 'pcg'); |
---|
85 | % |
---|
86 | |
---|
87 | % Author: Ashutosh Saxena <asaxena@cs.stanford.edu> |
---|
88 | % Modified by Min Sun <aliensun@stanford.edu> |
---|
89 | |
---|
90 | %DEBUG_FLAG = false; |
---|
91 | |
---|
92 | global A b S inq |
---|
93 | |
---|
94 | q = inq; |
---|
95 | |
---|
96 | % FEASIBLAE DATA CHECK |
---|
97 | [m,n] = size(A); % problem size: m examples, n features |
---|
98 | [l,n_S] = size(S); |
---|
99 | |
---|
100 | if ~isempty(S) |
---|
101 | if n_S ~= n; disp('size inconsistance'); return; end % stop if matrix dimension mismatch |
---|
102 | end |
---|
103 | |
---|
104 | %------------------------------------------------------------ |
---|
105 | % INITIALIZE |
---|
106 | %------------------------------------------------------------ |
---|
107 | |
---|
108 | % LOG BARRIER METHOD |
---|
109 | EPSILON_GAP = 5e-5; |
---|
110 | MU_t = 3;%15; %4; % for t -- log barrier. Changed ASH |
---|
111 | if(isempty(t_0)) t_0 = 500; end%500; end |
---|
112 | t = t_0; |
---|
113 | |
---|
114 | % SIGMOID APPROXIMATION |
---|
115 | NUMERICAL_LIMIT_EXP = 3e2; % this value is calculated such that exp(2*NUMERICAL_LIMIT_EXP) < inf |
---|
116 | MU_alpha = 1; ;%1.1; %1.5; % for sigmoidal |
---|
117 | MU_alpha2 = 2; %5;%1.5; %1.5; % for sigmoidal |
---|
118 | ALPHA_MAX = 3000; |
---|
119 | if(isempty(alpha_0 )) alpha_0 = 1e-1; end % ERROR, actually, it should be related to min(A*x-b) is set later in the code. See below |
---|
120 | ALPHA_MIN = 500.0; % should be atleast 100 -- ERROR |
---|
121 | |
---|
122 | % NEWTON PARAMETERS |
---|
123 | MAX_TNT_ITER = 50; %100; %25; % maximum Newton iteration |
---|
124 | %ABSTOL = 1e-8; % terminates when the norm of gradient < ABSTOL |
---|
125 | EPSILON = 1e-6; %1e-6; % terminate when lambdasqr_by_2 < EPSILON |
---|
126 | |
---|
127 | % LINE SEARCH PARAMETERS |
---|
128 | ALPHA_LineSearch= 0.01; % minimum fraction of decrease in norm(gradient) |
---|
129 | BETA = 0.5; % stepsize decrease factor |
---|
130 | MAX_LS_ITER = 25; % maximum backtracking line search iteration |
---|
131 | %MIN_STEP_SIZE = 1e-20; |
---|
132 | s0 = 1; |
---|
133 | |
---|
134 | |
---|
135 | |
---|
136 | % VARIABLE INITIALZE |
---|
137 | % find a feasible set for Sw+q <= 0 |
---|
138 | if isempty(Para.ptry) Para.ptry = 2:3:n; end % might not be true |
---|
139 | if isempty(Para.ptrz) Para.ptrz= 3:3:n; end% might not be true |
---|
140 | if isempty(Para.Dist_Start) Para.Dist_Start = (n/3+1); end |
---|
141 | w = zeros( n,1); |
---|
142 | w(Para.ptry) = -0.1; |
---|
143 | w(Para.ptrz) = 1; |
---|
144 | ll = S*w; |
---|
145 | factor = max(ll( Para.Dist_Start:end)); |
---|
146 | factor = 0.9/factor; |
---|
147 | w = w*factor; |
---|
148 | q_tilde = q + 1e-15; |
---|
149 | %q_solveInitial = q + 1e-10; |
---|
150 | % w = -S \ q_solveInitial; |
---|
151 | |
---|
152 | % if max( S*w+q_tilde) >= 0 |
---|
153 | % w0 = ones(n,1); |
---|
154 | % % find strickly feasible starting w |
---|
155 | % opt = sdpsettings('solver','sedumi','cachesolvers',1, 'verbose', 0); |
---|
156 | % w = sdpvar(n,1); |
---|
157 | % Strick_feasible_gap = (1/Para.ClosestDist -1/Para.FarestDist)/4*ones(size(q)); |
---|
158 | % Strick_feasible_gap(q == 0) = 1; |
---|
159 | % sol = solvesdp(set(S*w+q + Strick_feasible_gap<=0),norm(w0 - w),opt); |
---|
160 | % w = double(w); |
---|
161 | % end |
---|
162 | |
---|
163 | %save initial_w.mat w |
---|
164 | % load initial_w.mat |
---|
165 | |
---|
166 | % setting starting alpha_0 |
---|
167 | alpha_0 = 1 / max( abs(A*w-b) ); |
---|
168 | alfa = alpha_0; |
---|
169 | alfa_max = alfa; |
---|
170 | |
---|
171 | if max( S*w+q_tilde) >= 0 |
---|
172 | % disp('INFEASIBLE Start'); |
---|
173 | w0 = ones(n,1); |
---|
174 | % find strickly feasible starting w |
---|
175 | opt = sdpsettings('solver','sedumi','cachesolvers',1, 'verbose', 0); |
---|
176 | w = sdpvar(n,1); |
---|
177 | Strick_feasible_gap = (1/Para.ClosestDist -1/Para.FarestDist)/4*ones(size(q)); |
---|
178 | Strick_feasible_gap(q == 0) = 1; |
---|
179 | sol = solvesdp(set(S*w+q + Strick_feasible_gap<=0),norm(w0 - w),opt); |
---|
180 | w = double(w); |
---|
181 | end |
---|
182 | |
---|
183 | % ----------------------------------------------------------- |
---|
184 | % META LOOP FOR CHANGING t and ALPHA |
---|
185 | % ---------------------------------------------------------- |
---|
186 | history=[]; |
---|
187 | T_nt_hist=[]; |
---|
188 | |
---|
189 | if VERBOSE |
---|
190 | disp(sprintf('%15s %15s %11s %12s %10s %10s %7s %15s %15s',... |
---|
191 | 'exact obj', 'primal obj', 'alpha', 'MlogExpTerm',... |
---|
192 | 'norm(g)', 'lambda^2By2', ... |
---|
193 | 'gap', 'search_step', 'newton_steps')); |
---|
194 | %'normg_sigmoid', 'normg_t', ... |
---|
195 | end |
---|
196 | |
---|
197 | status = -1; |
---|
198 | %for outer_iter = 1:MAX_OUTER_ITER |
---|
199 | while (status~=2) |
---|
200 | |
---|
201 | history_Inner = []; |
---|
202 | |
---|
203 | dw = zeros(n,1); % dw newton step |
---|
204 | |
---|
205 | % Makes the profiler not work! |
---|
206 | % if (VERBOSE >= 2) |
---|
207 | % disp(sprintf('%s %15s %15s %11s %12s %10s %10s %10s %10s %9s %6s %6s %10s %10s %10s',... |
---|
208 | % 'iter', 'exact obj', 'primal obj', 'alpha', 'MlogExpTerm',... |
---|
209 | % 'In_step', 'stepsize','norm(g)', 'norm(dw)', ... |
---|
210 | % 'lambda^2By2','normg_sigmoid','normg_t', ... |
---|
211 | % 'normdw_t', 'normdw_sigmoid')); |
---|
212 | % end |
---|
213 | |
---|
214 | %------------------------------------------------------------ |
---|
215 | % MAIN LOOP |
---|
216 | %------------------------------------------------------------ |
---|
217 | total_linesearch = 0; |
---|
218 | ntiter = 0; |
---|
219 | %initialize Newton step |
---|
220 | logExpTerm = alfa*(A*w-b)'; |
---|
221 | expTerm = exp(logExpTerm); |
---|
222 | %expTermNeg = exp(-logExpTerm); |
---|
223 | expTerm_Rec = 1./(1+expTerm); |
---|
224 | inequalityTerm = (S*w+q)'; |
---|
225 | |
---|
226 | expTermNeg_Rec = 1./(1+exp(-logExpTerm)); |
---|
227 | g_sigmoid = (expTermNeg_Rec - expTerm_Rec) * t; |
---|
228 | gradphi_sigmoid = (g_sigmoid*A)'; %A1endPt) + ... |
---|
229 | g_t = (-1./inequalityTerm); % log barrier |
---|
230 | gradphi_t = (g_t*S)'; |
---|
231 | gradphi = (gradphi_sigmoid + gradphi_t); |
---|
232 | |
---|
233 | newtonLoop = (ntiter <= MAX_TNT_ITER); |
---|
234 | while newtonLoop |
---|
235 | ntiter = ntiter + 1; |
---|
236 | |
---|
237 | indicesInf_Pos = find(logExpTerm > NUMERICAL_LIMIT_EXP); |
---|
238 | indicesInf_Neg = find(logExpTerm < -NUMERICAL_LIMIT_EXP); |
---|
239 | %using logical mask instead to use 'find' might be faster. |
---|
240 | |
---|
241 | %indicesOkay = setdiff(1:length(logExpTerm), [indicesInf_Pos, indicesInf_Neg]); |
---|
242 | |
---|
243 | h_sigmoid = t * expTerm .* (expTerm_Rec.^2); |
---|
244 | h_sigmoid([indicesInf_Pos indicesInf_Neg]) = 0; |
---|
245 | %zeros(1, length(indicesInf_Pos)+length(indicesInf_Neg)); |
---|
246 | |
---|
247 | |
---|
248 | %------------------------------------------------------------ |
---|
249 | % CALCULATE NEWTON STEP |
---|
250 | %------------------------------------------------------------ |
---|
251 | |
---|
252 | hessphi_sigmoid = (sparse(1:m,1:m,h_sigmoid)*A)' * A; |
---|
253 | % hessphi_sigmoid = (h_sigmoid*(A.^2))'; %diagonal element only |
---|
254 | |
---|
255 | h_t = g_t.^2;% log barrier |
---|
256 | |
---|
257 | % if DEBUG_FLAG && (any(h_sigmoid<0) || any(h_t<0) ) |
---|
258 | % disp(sprintf('h negative = %5.4e; h_t negative = %5.4e',h_sigmoid(h_sigmoid<0), h_t(h_t<0))); |
---|
259 | % end |
---|
260 | |
---|
261 | |
---|
262 | hessphi_t = ( (sparse(1:length(h_t), 1:length(h_t), h_t) * S )' * S ) / (2*alfa); |
---|
263 | % hessphi_t = (h_t*(S.^2))' / (2*alfa); %Diagonal Element only |
---|
264 | |
---|
265 | %condest(hessphi_sigmoid) |
---|
266 | %condest(hessphi_t) |
---|
267 | hessphi = hessphi_sigmoid + hessphi_t; %HACK |
---|
268 | |
---|
269 | %condest(hessphi) |
---|
270 | dw = - hessphi \ gradphi; |
---|
271 | |
---|
272 | % dw = - (1./hessphi) .* gradphi; % Diagonal Hessian |
---|
273 | dw = dw / (2*alfa); |
---|
274 | |
---|
275 | % newton decrement=========================================== |
---|
276 | lambdasqr = full(-( gradphi'*dw) ); |
---|
277 | %lambdasqr2 = full(dw'*(hessphi+ hessphi_t)'*dw); |
---|
278 | % =========================================================== |
---|
279 | |
---|
280 | |
---|
281 | %------------------------------------------------------------ |
---|
282 | % BACKTRACKING LINE SEARCH |
---|
283 | %------------------------------------------------------------ |
---|
284 | s = s0; |
---|
285 | |
---|
286 | %maxInequality = max( S*(w+s*dw)+q_tilde); |
---|
287 | %while maxInequality >= 0 |
---|
288 | while any( ( S*(w+s*dw)+q_tilde) >= 0 ) |
---|
289 | s = BETA*s; |
---|
290 | %maxInequality = max( S*(w+s*dw)+q_tilde); |
---|
291 | end % first set the new w inside |
---|
292 | |
---|
293 | normg = norm(gradphi); |
---|
294 | lsiter = 0; |
---|
295 | backIterationLoop = true; |
---|
296 | while backIterationLoop |
---|
297 | lsiter = lsiter+1; |
---|
298 | new_w = w+s*dw; |
---|
299 | logExpTerm = alfa*(A*new_w-b)'; |
---|
300 | expTerm = exp(logExpTerm); |
---|
301 | |
---|
302 | % evaluate the gradphi |
---|
303 | inequalityTerm = (S*new_w+q)'; |
---|
304 | |
---|
305 | expTerm_Rec = 1./(1+expTerm); |
---|
306 | expTermNeg_Rec = 1./(1+exp(-logExpTerm)); |
---|
307 | g_sigmoid = (expTermNeg_Rec - expTerm_Rec) * t; |
---|
308 | gradphi_sigmoid = (g_sigmoid*A)'; %A1endPt) + ... |
---|
309 | g_t = (-1./inequalityTerm);% log barrier |
---|
310 | gradphi_t = (g_t*S)'; |
---|
311 | gradphi = (gradphi_sigmoid + gradphi_t); |
---|
312 | |
---|
313 | backIterationLoop = (lsiter <= MAX_LS_ITER) && ( norm(gradphi) > (1-ALPHA_LineSearch*s)*normg); |
---|
314 | s = BETA*s; |
---|
315 | end |
---|
316 | |
---|
317 | total_linesearch = total_linesearch + lsiter; |
---|
318 | |
---|
319 | % if VERBOSE >= 2 |
---|
320 | % Exact_phi_sigmoid = norm( logExpTerm/alfa, 1); % only for debug |
---|
321 | % Exact_f = Exact_phi_sigmoid + logBarrierValue / t; % MIN -- fixed |
---|
322 | % %normdw = norm(dw); |
---|
323 | % normg_sigmoid = norm(gradphi_sigmoid); |
---|
324 | % normg_t = norm(gradphi_t); |
---|
325 | % normdw_sigmoid = 0;%norm( dw_sigmoid); |
---|
326 | % normdw_t = 0;%norm( dw_t); |
---|
327 | % normdw = norm( dw); |
---|
328 | % f_new = phi_sigmoid + logBarrierValue / t; % MIN -- fixed ? |
---|
329 | % disp(sprintf('%4d %15.6e %15.6e %4.5e %10.2e %10.2e %10.2e %6d %6d %10.2e %6d %6d %6d %6d',... |
---|
330 | % ntiter,Exact_f, f_new, alfa, max(abs( logExpTerm)), s0, s, normg, normdw, lambdasqr/2, ... |
---|
331 | % normg_sigmoid, normg_t, normdw_t, normdw_sigmoid)); |
---|
332 | % end |
---|
333 | |
---|
334 | %------------------------------------------------------------ |
---|
335 | % STOPPING CRITERION |
---|
336 | %------------------------------------------------------------ |
---|
337 | |
---|
338 | if (lambdasqr/2 <= EPSILON) |
---|
339 | status = 1; |
---|
340 | %disp('Minimal Newton decrement reached'); |
---|
341 | end |
---|
342 | |
---|
343 | newtonLoop = (ntiter <= MAX_TNT_ITER) && (lambdasqr/2 > EPSILON) && (lsiter < MAX_LS_ITER); % |
---|
344 | |
---|
345 | % set new w |
---|
346 | w = new_w; |
---|
347 | %alfa = min(alfa*MU_alpha, alfa_max); |
---|
348 | % If you want to do the above, i.e., increase the value in each |
---|
349 | % Newton step, uncomment the gradient calculation, if ntiter == 1 |
---|
350 | |
---|
351 | end % -----------END of the MAIN LOOP----------------------------- |
---|
352 | |
---|
353 | % Tighten the sigmoid and log approximation |
---|
354 | |
---|
355 | |
---|
356 | gap = m/t; |
---|
357 | if VERBOSE |
---|
358 | Exact_phi_sigmoid = norm( logExpTerm/alfa, 1); % only for debug |
---|
359 | |
---|
360 | normg_sigmoid = norm(gradphi_sigmoid); |
---|
361 | normg_t = norm(gradphi_t); |
---|
362 | normdw = norm( dw); |
---|
363 | normg = norm( gradphi); |
---|
364 | expTermNeg = exp(-logExpTerm); |
---|
365 | phi_sigmoid = ( log( 1+ expTermNeg ) + log( 1+ expTerm) ); |
---|
366 | phi_sigmoid( indicesInf_Pos ) = logExpTerm( indicesInf_Pos ); |
---|
367 | phi_sigmoid( indicesInf_Neg ) = -logExpTerm( indicesInf_Neg ); |
---|
368 | phi_sigmoid = (1/alfa) * sum(phi_sigmoid,2); |
---|
369 | |
---|
370 | logBarrierValue = -sum(log( -inequalityTerm ) ); |
---|
371 | |
---|
372 | f_new = phi_sigmoid + logBarrierValue / t; % MIN -- fixed ? |
---|
373 | Exact_f = Exact_phi_sigmoid + logBarrierValue / t; % MIN -- fixed |
---|
374 | disp(sprintf('%15.6e %15.6e %4.5e %10.2e %10.2e %10.2e %10.2d %10d %10d',... |
---|
375 | Exact_f, f_new, alfa, max(abs( logExpTerm)), normg, lambdasqr/2, ... |
---|
376 | gap, total_linesearch, ntiter)); |
---|
377 | % normg_sigmoid, normg_t, ... |
---|
378 | end |
---|
379 | |
---|
380 | T_nt_hist = [T_nt_hist history_Inner]; |
---|
381 | history=[history [length(history_Inner); gap]]; |
---|
382 | |
---|
383 | %if alfa >= ALPHA_MIN/MU_alpha |
---|
384 | epsilon_gap = EPSILON_GAP; |
---|
385 | %else |
---|
386 | % epsilon_gap = (1-alfa/ALPHA_MIN)*1e2 + (alfa/ALPHA_MIN)*EPSILON_GAP; |
---|
387 | %end |
---|
388 | |
---|
389 | if (alfa > ALPHA_MIN) && (gap < epsilon_gap) && (status >= 1) |
---|
390 | %disp('Alpha and Gap reached'); |
---|
391 | status=2; |
---|
392 | end |
---|
393 | % alfa = min(alfa*MU_alpha, ALPHA_MAX); |
---|
394 | % alfa_max = min(alfa_max*MU_alpha2, ALPHA_MAX); |
---|
395 | % if (alfa > ALPHA_MIN) && (status >= 1) |
---|
396 | % t = MU_t*t; % if just waiting for gap, then increase it at double speed |
---|
397 | % end |
---|
398 | |
---|
399 | t = MU_t*t; |
---|
400 | alfa = min(alfa*MU_alpha2, ALPHA_MAX); |
---|
401 | alfa_max = min(alfa*MU_alpha2, ALPHA_MAX); |
---|
402 | |
---|
403 | end |
---|
404 | |
---|
405 | return; |
---|
406 | |
---|