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, A, b, S, q, method, ptol, pmaxi, VERBOSE) |
---|
41 | % |
---|
42 | % l1 penalty approximate Problem with inequality constrain Solver |
---|
43 | % |
---|
44 | % exact problems form: |
---|
45 | % |
---|
46 | % minimize norm( Aw - b, 1) st, Sw+q<=0 |
---|
47 | % |
---|
48 | % sigmoid approximate form: |
---|
49 | % |
---|
50 | % minimize sum_i (1/alfa)*( log( 1+exp(-alfa*( A(i,:)*w - b(i)))) + ... |
---|
51 | % log( 1+exp(alfa*( A(i,:)*w - b(i)))) ) |
---|
52 | % st, Sw+q <= 0 (log barrier approximation) |
---|
53 | % |
---|
54 | % where variable is w and problem data are A, b, S, and q. |
---|
55 | % |
---|
56 | % INPUT |
---|
57 | % |
---|
58 | % A : mxn matrix; |
---|
59 | % b : m vector; |
---|
60 | % S : lxn matrix; |
---|
61 | % q : l vector; |
---|
62 | % |
---|
63 | % method : string; search direction method type |
---|
64 | % 'cg' : conjugate gradients method, 'pcg' |
---|
65 | % 'pcg' : preconditioned conjugate gradients method |
---|
66 | % 'exact': exact method (default value) |
---|
67 | % ptol : scalar; pcg relative tolerance. if empty, use adaptive rule. |
---|
68 | % pmaxi : scalar: pcg maximum iteration. if empty, use default value (500).a |
---|
69 | % mu : approximate_closeness factor; bigger -> more accurate |
---|
70 | % VERBOSE : enable disp |
---|
71 | % |
---|
72 | % OUTPUT |
---|
73 | % |
---|
74 | % w : n vector; classifier |
---|
75 | % status : scalar; +1: success, -1: maxiter exceeded |
---|
76 | % history : |
---|
77 | % row 1) phi (objective value) |
---|
78 | % row 2) norm(gradient of phi) |
---|
79 | % row 3) cumulative cg iterations |
---|
80 | % |
---|
81 | % USAGE EXAMPLE |
---|
82 | % |
---|
83 | % [w,status] = (A, b, S, q, 'pcg'); |
---|
84 | % |
---|
85 | |
---|
86 | % template by Kwangmoo Koh <deneb1@stanford.edu> |
---|
87 | % Modified by Min Sun <aliensun@stanford.edu> |
---|
88 | |
---|
89 | %DEBUG_FLAG = false; |
---|
90 | S_temp = S; |
---|
91 | A_temp = A; |
---|
92 | global A S h_tt h_sigmoidt A2 S2; |
---|
93 | S = S_temp; |
---|
94 | A = A_temp; |
---|
95 | clear A_temp S_temp |
---|
96 | A2 = (A.^2)'; |
---|
97 | S2 = (S.^2)'; |
---|
98 | |
---|
99 | % FEASIBLAE DATA CHECK |
---|
100 | [m,n] = size(A); % problem size: m examples, n features |
---|
101 | [l,n_S] = size(S); |
---|
102 | |
---|
103 | if ~isempty(S) |
---|
104 | if n_S ~= n; disp('size inconsistance'); return; end % stop if matrix dimension mismatch |
---|
105 | end |
---|
106 | |
---|
107 | |
---|
108 | %------------------------------------------------------------ |
---|
109 | % INITIALIZE |
---|
110 | %------------------------------------------------------------ |
---|
111 | |
---|
112 | % LOG BARRIER METHOD |
---|
113 | EPSILON_GAP = 5e-5; |
---|
114 | MU_t = 3;%15; %4; % for t -- log barrier. Changed ASH |
---|
115 | if(isempty(t_0)) t_0 = 1; end%500; end |
---|
116 | t = t_0; |
---|
117 | |
---|
118 | % SIGMOID APPROXIMATION |
---|
119 | NUMERICAL_LIMIT_EXP = 3e2; % this value is calculated such that exp(2*NUMERICAL_LIMIT_EXP) < inf |
---|
120 | MU_alpha = 1; ;%1.1; %1.5; % for sigmoidal |
---|
121 | MU_alpha2 = 3; %5;%1.5; %1.5; % for sigmoidal |
---|
122 | ALPHA_MAX = 500; |
---|
123 | 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 |
---|
124 | ALPHA_MIN = 200.0; % should be atleast 100 -- ERROR |
---|
125 | |
---|
126 | % NEWTON PARAMETERS |
---|
127 | MAX_TNT_ITER = 50; %100; %25; % maximum Newton iteration |
---|
128 | %ABSTOL = 1e-8; % terminates when the norm of gradient < ABSTOL |
---|
129 | EPSILON = 1e-6; %1e-6; % terminate when lambdasqr_by_2 < EPSILON |
---|
130 | |
---|
131 | % LINE SEARCH PARAMETERS |
---|
132 | ALPHA_LineSearch= 0.01; % minimum fraction of decrease in norm(gradient) |
---|
133 | BETA = 0.5; % stepsize decrease factor |
---|
134 | MAX_LS_ITER = 100; % maximum backtracking line search iteration |
---|
135 | %MIN_STEP_SIZE = 1e-20; |
---|
136 | s0 = 1; |
---|
137 | |
---|
138 | % PCG parameters |
---|
139 | % PCG METHOD |
---|
140 | if(isempty(pmaxi)) pcgmaxi = 100; else pcgmaxi = pmaxi; end |
---|
141 | if(isempty(ptol )) pcgtol = 1e-4; else pcgtol = ptol; end |
---|
142 | |
---|
143 | |
---|
144 | |
---|
145 | % VARIABLE INITIALZE |
---|
146 | % find a feasible set for Sw+q <= 0 |
---|
147 | q_tilde = q + 1e-15; |
---|
148 | q_solveInitial = q + 1e-10; |
---|
149 | w = -S \ q_solveInitial; |
---|
150 | if isempty(w) |
---|
151 | w0 = ones(n,1); |
---|
152 | % find strickly feasible starting w |
---|
153 | opt = sdpsettings('solver','sedumi','cachesolvers',1, 'verbose', 0); |
---|
154 | w = sdpvar(n,1); |
---|
155 | Strick_feasible_gap = (1/Para.ClosestDist -1/Para.FarestDist)/4*ones(size(q)); |
---|
156 | Strick_feasible_gap(q == 0) = 1; |
---|
157 | sol = solvesdp(set(S*w+q + Strick_feasible_gap<=0),norm(w0 - w),opt); |
---|
158 | w = double(w); |
---|
159 | end |
---|
160 | |
---|
161 | % setting starting alpha_0 |
---|
162 | alpha_0 = 1 / max( abs(A*w-b) ); |
---|
163 | alfa = alpha_0; |
---|
164 | alfa_max = alfa; |
---|
165 | |
---|
166 | if max( S*w+q_tilde) >= 0 |
---|
167 | disp('INFEASIBLE Start'); |
---|
168 | end |
---|
169 | |
---|
170 | % ----------------------------------------------------------- |
---|
171 | % META LOOP FOR CHANGING t and ALPHA |
---|
172 | % ---------------------------------------------------------- |
---|
173 | history=[]; |
---|
174 | T_nt_hist=[]; |
---|
175 | |
---|
176 | if VERBOSE |
---|
177 | disp(sprintf('%15s %15s %11s %12s %10s %10s %10s %10s %10s %7s %15s %15s',... |
---|
178 | 'exact obj', 'primal obj', 'alpha', 'MlogExpTerm',... |
---|
179 | 'norm(g)', 'norm(dw)', 'lambda^2By2', 'normg_sigmoid', 'normg_t', 'gap', ... |
---|
180 | 'search_step', 'newton_steps')); |
---|
181 | end |
---|
182 | |
---|
183 | status = -1; |
---|
184 | %for outer_iter = 1:MAX_OUTER_ITER |
---|
185 | while (status~=2) |
---|
186 | |
---|
187 | history_Inner = []; |
---|
188 | |
---|
189 | dw = zeros(n,1); % dw newton step |
---|
190 | |
---|
191 | % Makes the profiler not work! |
---|
192 | % if (VERBOSE >= 2) |
---|
193 | % disp(sprintf('%s %15s %15s %11s %12s %10s %10s %10s %10s %9s %6s %6s %10s %10s %10s',... |
---|
194 | % 'iter', 'exact obj', 'primal obj', 'alpha', 'MlogExpTerm',... |
---|
195 | % 'In_step', 'stepsize','norm(g)', 'norm(dw)', ... |
---|
196 | % 'lambda^2By2','normg_sigmoid','normg_t', ... |
---|
197 | % 'normdw_t', 'normdw_sigmoid')); |
---|
198 | % end |
---|
199 | |
---|
200 | %------------------------------------------------------------ |
---|
201 | % MAIN LOOP |
---|
202 | %------------------------------------------------------------ |
---|
203 | total_linesearch = 0; |
---|
204 | ntiter = 0; |
---|
205 | %initialize Newton step |
---|
206 | logExpTerm = alfa*(A*w-b)'; |
---|
207 | expTerm = exp(logExpTerm); |
---|
208 | %expTermNeg = exp(-logExpTerm); |
---|
209 | expTerm_Rec = 1./(1+expTerm); |
---|
210 | inequalityTerm = (S*w+q)'; |
---|
211 | |
---|
212 | expTermNeg_Rec = 1./(1+exp(-logExpTerm)); |
---|
213 | g_sigmoid = (expTermNeg_Rec - expTerm_Rec) * t; |
---|
214 | gradphi_sigmoid = (g_sigmoid*A)'; %A1endPt) + ... |
---|
215 | g_t = sparse(-1./inequalityTerm); % log barrier |
---|
216 | gradphi_t = (g_t*S)'; |
---|
217 | gradphi = (gradphi_sigmoid + gradphi_t); |
---|
218 | |
---|
219 | newtonLoop = (ntiter <= MAX_TNT_ITER); |
---|
220 | while newtonLoop |
---|
221 | ntiter = ntiter + 1; |
---|
222 | |
---|
223 | indicesInf_Pos = find(logExpTerm > NUMERICAL_LIMIT_EXP); |
---|
224 | indicesInf_Neg = find(logExpTerm < -NUMERICAL_LIMIT_EXP); |
---|
225 | %indicesOkay = setdiff(1:length(logExpTerm), [indicesInf_Pos, indicesInf_Neg]); |
---|
226 | |
---|
227 | h_sigmoid = t * expTerm .* (expTerm_Rec.^2); |
---|
228 | h_sigmoid([indicesInf_Pos indicesInf_Neg]) = zeros(1, length(indicesInf_Pos)+length(indicesInf_Neg)); |
---|
229 | |
---|
230 | h_t = g_t.^2;% log barrier |
---|
231 | |
---|
232 | norm(gradphi) |
---|
233 | %------------------------------------------------------------ |
---|
234 | % CALCULATE NEWTON STEP |
---|
235 | %------------------------------------------------------------ |
---|
236 | switch lower(method) |
---|
237 | case 'pcg' |
---|
238 | %M_cond = sparse(1:size(A2,1), 1:size(A2,1), 1./(A2*h_sigmoid'+S2*h_t') ); |
---|
239 | h_sigmoidt = h_sigmoid'; |
---|
240 | h_tt = h_t'; |
---|
241 | if (isempty(ptol)) pcgtol = min(0.001,norm(gradphi)); end |
---|
242 | [dw, pflg, prelres, pitr, presvec] = ... |
---|
243 | pcg(@AXfunc,-gradphi,pcgtol,pcgmaxi, @Mfunc); |
---|
244 | % pcg(@AXfunc,-gradphi,pcgtol,pcgmaxi,@Mfunc,[],[],... |
---|
245 | % A, h_sigmoid, S, h_t, 1./(A2'*h_sigmoid+S2'*h_t)); |
---|
246 | if (pitr == 0) pitr = pcgmaxi; end |
---|
247 | otherwise |
---|
248 | hessphi_sigmoid = (sparse(1:m,1:m,h_sigmoid)*A)' * A; |
---|
249 | |
---|
250 | hessphi_t = ( (sparse(1:length(h_t), 1:length(h_t), h_t) * S )' * S ) / (2*alfa); |
---|
251 | |
---|
252 | hessphi = hessphi_sigmoid + hessphi_t; %HACK |
---|
253 | |
---|
254 | dw = - hessphi \ gradphi; |
---|
255 | end |
---|
256 | % dw = - (1./hessphi) .* gradphi; % Diagonal Hessian |
---|
257 | dw = dw / (2*alfa); |
---|
258 | |
---|
259 | % newton decrement=========================================== |
---|
260 | lambdasqr = full(-( gradphi'*dw) ); |
---|
261 | %lambdasqr2 = full(dw'*(hessphi+ hessphi_t)'*dw); |
---|
262 | % =========================================================== |
---|
263 | |
---|
264 | |
---|
265 | %------------------------------------------------------------ |
---|
266 | % BACKTRACKING LINE SEARCH |
---|
267 | %------------------------------------------------------------ |
---|
268 | s = s0; |
---|
269 | |
---|
270 | maxInequality = max( S*(w+s*dw)+q_tilde); |
---|
271 | while maxInequality >= 0 |
---|
272 | maxInequality = max( S*(w+s*dw)+q_tilde); |
---|
273 | s = BETA*s; |
---|
274 | end % first set the new w inside |
---|
275 | |
---|
276 | normg = norm(gradphi); |
---|
277 | lsiter = 0; |
---|
278 | backIterationLoop = true; |
---|
279 | while backIterationLoop |
---|
280 | lsiter = lsiter+1; |
---|
281 | new_w = w+s*dw; |
---|
282 | logExpTerm = alfa*(A*new_w-b)'; |
---|
283 | expTerm = exp(logExpTerm); |
---|
284 | |
---|
285 | % evaluate the gradphi |
---|
286 | inequalityTerm = (S*new_w+q)'; |
---|
287 | |
---|
288 | expTerm_Rec = 1./(1+expTerm); |
---|
289 | expTermNeg_Rec = 1./(1+exp(-logExpTerm)); |
---|
290 | g_sigmoid = (expTermNeg_Rec - expTerm_Rec) * t; |
---|
291 | gradphi_sigmoid = (g_sigmoid*A)'; %A1endPt) + ... |
---|
292 | g_t = sparse(-1./inequalityTerm);% log barrier |
---|
293 | gradphi_t = (g_t*S)'; |
---|
294 | gradphi = (gradphi_sigmoid + gradphi_t); |
---|
295 | |
---|
296 | backIterationLoop = (lsiter <= MAX_LS_ITER) && ( norm(gradphi) > (1-ALPHA_LineSearch*s)*normg); |
---|
297 | s = BETA*s; |
---|
298 | end |
---|
299 | |
---|
300 | total_linesearch = total_linesearch + lsiter; |
---|
301 | |
---|
302 | % if VERBOSE >= 2 |
---|
303 | % Exact_phi_sigmoid = norm( logExpTerm/alfa, 1); % only for debug |
---|
304 | % Exact_f = Exact_phi_sigmoid + logBarrierValue / t; % MIN -- fixed |
---|
305 | % %normdw = norm(dw); |
---|
306 | % normg_sigmoid = norm(gradphi_sigmoid); |
---|
307 | % normg_t = norm(gradphi_t); |
---|
308 | % normdw_sigmoid = 0;%norm( dw_sigmoid); |
---|
309 | % normdw_t = 0;%norm( dw_t); |
---|
310 | % normdw = norm( dw); |
---|
311 | % f_new = phi_sigmoid + logBarrierValue / t; % MIN -- fixed ? |
---|
312 | % disp(sprintf('%4d %15.6e %15.6e %4.5e %10.2e %10.2e %10.2e %6d %6d %10.2e %6d %6d %6d %6d',... |
---|
313 | % ntiter,Exact_f, f_new, alfa, max(abs( logExpTerm)), s0, s, normg, normdw, lambdasqr/2, ... |
---|
314 | % normg_sigmoid, normg_t, normdw_t, normdw_sigmoid)); |
---|
315 | % end |
---|
316 | |
---|
317 | %------------------------------------------------------------ |
---|
318 | % STOPPING CRITERION |
---|
319 | %------------------------------------------------------------ |
---|
320 | |
---|
321 | if (lambdasqr/2 <= EPSILON) |
---|
322 | status = 1; |
---|
323 | %disp('Minimal Newton decrement reached'); |
---|
324 | end |
---|
325 | |
---|
326 | lsiter |
---|
327 | MAX_LS_ITER |
---|
328 | newtonLoop = (ntiter <= MAX_TNT_ITER) && (lambdasqr/2 > EPSILON) && (lsiter < MAX_LS_ITER); % |
---|
329 | |
---|
330 | % set new w |
---|
331 | w = new_w; |
---|
332 | %alfa = min(alfa*MU_alpha, alfa_max); |
---|
333 | % If you want to do the above, i.e., increase the value in each |
---|
334 | % Newton step, uncomment the gradient calculation, if ntiter == 1 |
---|
335 | |
---|
336 | end % -----------END of the MAIN LOOP----------------------------- |
---|
337 | |
---|
338 | % Tighten the sigmoid and log approximation |
---|
339 | |
---|
340 | |
---|
341 | gap = m/t; |
---|
342 | if VERBOSE |
---|
343 | Exact_phi_sigmoid = norm( logExpTerm/alfa, 1); % only for debug |
---|
344 | |
---|
345 | normg_sigmoid = norm(gradphi_sigmoid); |
---|
346 | normg_t = norm(gradphi_t); |
---|
347 | normdw = norm( dw); |
---|
348 | normg = norm( gradphi); |
---|
349 | expTermNeg = exp(-logExpTerm); |
---|
350 | phi_sigmoid = ( log( 1+ expTermNeg ) + log( 1+ expTerm) ); |
---|
351 | phi_sigmoid( indicesInf_Pos ) = logExpTerm( indicesInf_Pos ); |
---|
352 | phi_sigmoid( indicesInf_Neg ) = -logExpTerm( indicesInf_Neg ); |
---|
353 | phi_sigmoid = (1/alfa) * sum(phi_sigmoid,2); |
---|
354 | |
---|
355 | logBarrierValue = -sum(log( -inequalityTerm ) ); |
---|
356 | |
---|
357 | f_new = phi_sigmoid + logBarrierValue / t; % MIN -- fixed ? |
---|
358 | Exact_f = Exact_phi_sigmoid + logBarrierValue / t; % MIN -- fixed |
---|
359 | disp(sprintf('%15.6e %15.6e %4.5e %10.2e %10.2e %12.2e %10.2e %6d %6d %10.2d %10d %10d',... |
---|
360 | Exact_f, f_new, alfa, max(abs( logExpTerm)), normg, normdw, lambdasqr/2, ... |
---|
361 | normg_sigmoid, normg_t, gap, total_linesearch, ntiter)); |
---|
362 | end |
---|
363 | |
---|
364 | T_nt_hist = [T_nt_hist history_Inner]; |
---|
365 | history=[history [length(history_Inner); gap]]; |
---|
366 | |
---|
367 | %if alfa >= ALPHA_MIN/MU_alpha |
---|
368 | epsilon_gap = EPSILON_GAP; |
---|
369 | %else |
---|
370 | % epsilon_gap = (1-alfa/ALPHA_MIN)*1e2 + (alfa/ALPHA_MIN)*EPSILON_GAP; |
---|
371 | %end |
---|
372 | |
---|
373 | if (alfa > ALPHA_MIN) && (gap < epsilon_gap) && (status >= 1) |
---|
374 | disp('Alpha and Gap reached'); |
---|
375 | status=2; |
---|
376 | end |
---|
377 | % alfa = min(alfa*MU_alpha, ALPHA_MAX); |
---|
378 | % alfa_max = min(alfa_max*MU_alpha2, ALPHA_MAX); |
---|
379 | if (alfa > ALPHA_MIN) && (status >= 1) |
---|
380 | t = MU_t*t; % if just waiting for gap, then increase it at double speed |
---|
381 | end |
---|
382 | |
---|
383 | t = MU_t*t; |
---|
384 | alfa = min(alfa*MU_alpha2, ALPHA_MAX); |
---|
385 | alfa_max = min(alfa*MU_alpha2, ALPHA_MAX); |
---|
386 | |
---|
387 | end |
---|
388 | |
---|
389 | return; |
---|
390 | |
---|
391 | |
---|
392 | |
---|
393 | %------------------------------------------------------------ |
---|
394 | % CALL BACK FUNCTIONS FOR PCG |
---|
395 | %------------------------------------------------------------ |
---|
396 | function y = AXfunc(x) |
---|
397 | global A S h_tt h_sigmoidt; |
---|
398 | % global hessphi; |
---|
399 | % y = hessphi*x; |
---|
400 | y = ((h_sigmoidt.*(A*x))'*A)'+... |
---|
401 | ((h_tt.*(S*x))'*S)'; |
---|
402 | return; |
---|
403 | |
---|
404 | function y = Mfunc(x) |
---|
405 | global h_tt h_sigmoidt A2 S2; |
---|
406 | y = x./(A2*h_sigmoidt+S2*h_tt); |
---|
407 | return; |
---|