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, alpha, status, history] = Min_minimizer( Para, t, alpha, w, A, b, S, q, method, ptol, pmaxi, mu, VERBOSE, LogBarrierFlg) |
---|
40 | % |
---|
41 | % l1 penalty approximate Problem with inequality constrain Solver |
---|
42 | % |
---|
43 | % exact problems form: |
---|
44 | % |
---|
45 | % minimize norm( Aw - b, 1) st, Sw<=q |
---|
46 | % |
---|
47 | % sigmoid approximate form: |
---|
48 | % |
---|
49 | % minimize sum_i (1/alpha)*( log( 1+exp(-alpha*( A(i,:)*w - b(i)))) + ... |
---|
50 | % log( 1+exp(alpha*( A(i,:)*w - b(i)))) ) |
---|
51 | % st, Sw <= q (log barrier approximation) |
---|
52 | % |
---|
53 | % where variable is w and problem data are A, b, S, and q. |
---|
54 | % |
---|
55 | % INPUT |
---|
56 | % |
---|
57 | % A : mxn matrix; |
---|
58 | % b : m vector; |
---|
59 | % S : lxn matrix; |
---|
60 | % q : l vector; |
---|
61 | % |
---|
62 | % method : string; search direction method type |
---|
63 | % 'cg' : conjugate gradients method, 'pcg' |
---|
64 | % 'pcg' : preconditioned conjugate gradients method |
---|
65 | % 'exact': exact method (default value) |
---|
66 | % ptol : scalar; pcg relative tolerance. if empty, use adaptive rule. |
---|
67 | % pmaxi : scalar: pcg maximum iteration. if empty, use default value (500).a |
---|
68 | % mu : approximate_closeness factor; bigger -> more accurate |
---|
69 | % VERBOSE : enable disp |
---|
70 | % |
---|
71 | % OUTPUT |
---|
72 | % |
---|
73 | % w : n vector; classifier |
---|
74 | % status : scalar; +1: success, -1: maxiter exceeded |
---|
75 | % history : |
---|
76 | % row 1) phi (objective value) |
---|
77 | % row 2) norm(gradient of phi) |
---|
78 | % row 3) cumulative cg iterations |
---|
79 | % |
---|
80 | % USAGE EXAMPLE |
---|
81 | % |
---|
82 | % [w,status] = (A, b, S, q, 'pcg'); |
---|
83 | % |
---|
84 | |
---|
85 | % template by Kwangmoo Koh <deneb1@stanford.edu> |
---|
86 | % Modified by Min Sun <aliensun@stanford.edu> |
---|
87 | %------------------------------------------------------------ |
---|
88 | % INITIALIZE |
---|
89 | %------------------------------------------------------------ |
---|
90 | |
---|
91 | % NEWTON PARAMETERS |
---|
92 | MAX_TNT_ITER = 200; % maximum (truncated) Newton iteration |
---|
93 | ABSTOL = 1e-8; % terminates when the norm of gradient < ABSTOL |
---|
94 | EPSILON = 1e0; %1e-6; % terminate when lambdasqr_by_2 < EPSILON |
---|
95 | |
---|
96 | % LINE SEARCH PARAMETERS |
---|
97 | ALPHA_LineSearch= 0.01; % minimum fraction of decrease in norm(gradient) |
---|
98 | BETA = 0.5; % stepsize decrease factor |
---|
99 | MAX_LS_ITER = 100; % maximum backtracking line search iteration |
---|
100 | |
---|
101 | if(isempty(pmaxi)) pcgmaxi = 500; else pcgmaxi = pmaxi; end |
---|
102 | if(isempty(ptol )) pcgtol = 1e-4; else pcgtol = ptol; end |
---|
103 | |
---|
104 | % sigmoid approximation |
---|
105 | ALPHA_MAX = 500000; |
---|
106 | Eps = 0; |
---|
107 | if(isempty(alpha )) alpha = 1; end |
---|
108 | %MIN_ALPHA = 10000; |
---|
109 | mu = 1; |
---|
110 | %if(isempty(mu)) mu = 1; end % how to tune mu so that norm(gradient) decrease??? |
---|
111 | |
---|
112 | % log barrier moethd |
---|
113 | h_t_thre = 1e15; |
---|
114 | %t = 1; |
---|
115 | %mu_t = 1.5; |
---|
116 | |
---|
117 | % Data related Parameter |
---|
118 | [m,n] = size(A); % problem size: m examples, n features |
---|
119 | [l,n_S] = size(S); |
---|
120 | |
---|
121 | if ~isempty(S) |
---|
122 | if n_S ~= n; disp('size inconsistance'); return; end % stop if matrix dimension mismatch |
---|
123 | end |
---|
124 | |
---|
125 | % INITIALIZE |
---|
126 | pobj = Inf; s = Inf; pitr = 0 ; pflg = 0 ; prelres = 0; pcgiter = 0; |
---|
127 | history = []; |
---|
128 | |
---|
129 | %w0 = ones(n,1); |
---|
130 | % find strickly feasible starting w |
---|
131 | %opt = sdpsettings('solver','sedumi','cachesolvers',1, 'verbose', 0); |
---|
132 | %w = sdpvar(n,1); |
---|
133 | %Strick_feasible_gap = (1/Para.ClosestDist -1/Para.FarestDist)/4*ones(size(q)); |
---|
134 | %Strick_feasible_gap(q == 0) = 1; |
---|
135 | %sol = solvesdp(set(S*w+q + Strick_feasible_gap<=0),norm(w0 - w),opt); |
---|
136 | %w = double(w); |
---|
137 | dw = zeros(n,1); % dw newton step |
---|
138 | |
---|
139 | if VERBOSE |
---|
140 | disp(sprintf('%s %15s %15s %11s %12s %6s %10s %6s %6s %6s %10s',... |
---|
141 | 'iter', 'exact obj', 'primal obj', 'alpha', 'MlogExpTerm', 'stepsize','norm(g)', 'lambda^2By2','p_flg','p_res','p_itr')); |
---|
142 | end |
---|
143 | |
---|
144 | %------------------------------------------------------------ |
---|
145 | % MAIN LOOP |
---|
146 | %------------------------------------------------------------ |
---|
147 | |
---|
148 | for ntiter = 0:MAX_TNT_ITER |
---|
149 | |
---|
150 | % gradient related |
---|
151 | logExpTerm = alpha*(A*w-b); |
---|
152 | expTerm = exp(logExpTerm); |
---|
153 | expTermNeg = exp(-logExpTerm); |
---|
154 | g = sparse( 1./(1+expTermNeg) - 1./(1+expTerm) )*t; |
---|
155 | if LogBarrierFlg |
---|
156 | g_t = sparse(-1./(S*w+q));% log barrier |
---|
157 | end |
---|
158 | |
---|
159 | % ======= origin of bad condition of Hessian |
---|
160 | h = exp( logExpTerm - 2*log(1+expTerm) )*t; |
---|
161 | h_cond = sum( h(1:Para.A1endPt) <=Eps ) >0 &&... |
---|
162 | sum( h( (Para.A1endPt+1):(Para.A2endPt)) <=Eps) >0 &&... |
---|
163 | sum( h( (Para.A2endPt+1):(Para.A3endPt)) <=Eps) >0 &&... |
---|
164 | sum( h( (Para.A3endPt+1):end) <=Eps) >0; |
---|
165 | if h_cond % heuristic of condition of hessphi (usually stop 1 to 5 step earlier) |
---|
166 | if VERBOSE |
---|
167 | % disp('bad condition'); |
---|
168 | end |
---|
169 | % break; |
---|
170 | end |
---|
171 | % ========================================== |
---|
172 | if LogBarrierFlg |
---|
173 | h_t = 1./(S*w+q).^2;% log barrier |
---|
174 | end |
---|
175 | if any(h_t >= h_t_thre) % heuristic of condition of hessphi (usually stop 1 to 5 step earlier) |
---|
176 | if VERBOSE |
---|
177 | % disp('g_t bad condition'); |
---|
178 | end |
---|
179 | % break; |
---|
180 | end |
---|
181 | |
---|
182 | % gradphi = A'*g; |
---|
183 | gradphi = A( 1:Para.A1endPt,:)'*g(1:Para.A1endPt) + ... |
---|
184 | A( (1+Para.A1endPt):Para.A2endPt,:)'*g(( Para.A1endPt+1):Para.A2endPt) + ... |
---|
185 | A( (1+Para.A2endPt):Para.A3endPt,:)'*g( (Para.A2endPt+1):Para.A3endPt) + ... |
---|
186 | A( (1+Para.A3endPt):end,:)'*g( (Para.A3endPt+1):end); |
---|
187 | if LogBarrierFlg |
---|
188 | gradphi_t = S'* g_t; |
---|
189 | end |
---|
190 | normg = norm(gradphi+gradphi_t); |
---|
191 | |
---|
192 | phi = sum( (1/alpha)*( log( 1+ expTermNeg) + log( 1+ expTerm)), 1); |
---|
193 | Exact_phi = norm( logExpTerm/alpha, 1); % only for debug |
---|
194 | logBarrierValue = -sum(log( -S*w-q)); |
---|
195 | %------------------------------------------------------------ |
---|
196 | % CALCULATE NEWTON STEP |
---|
197 | %------------------------------------------------------------ |
---|
198 | switch lower(method) |
---|
199 | case 'pcg' |
---|
200 | % if (isempty(ptol)) pcgtol = min(0.1,norm(gradphi)); end |
---|
201 | % [dw, pflg, prelres, pitr, presvec] = ... |
---|
202 | % pcg(@AXfunc,-gradphi,pcgtol,pcgmaxi,@Mfunc,[],[],... |
---|
203 | % A,h,2*lambda,1./(A2'*h+2*lambda)); |
---|
204 | % if (pitr == 0) pitr = pcgmaxi; end |
---|
205 | |
---|
206 | case 'cg' |
---|
207 | % if (isempty(ptol)) pcgtol = min(0.1,norm(gradphi)); end |
---|
208 | % [dw, pflg, prelres, pitr, presvec] = ... |
---|
209 | % pcg(@AXfunc,-gradphi,pcgtol,pcgmaxi,[],[],[],... |
---|
210 | % A,h,2*lambda,[]); |
---|
211 | % if (pitr == 0) pitr = pcgmaxi; end |
---|
212 | |
---|
213 | otherwise % exact method |
---|
214 | % hessphi = A'*sparse(1:m,1:m,h)*A; |
---|
215 | hessphi = ... |
---|
216 | ( A( 1:Para.A1endPt,:)' *... |
---|
217 | sparse(1:Para.A1endPt, 1:Para.A1endPt, h(1:Para.A1endPt)) * ... |
---|
218 | A( 1:Para.A1endPt,:) + ... |
---|
219 | A( (1+Para.A1endPt):Para.A2endPt,:)' *... |
---|
220 | sparse(1:(Para.A2endPt-Para.A1endPt), 1:(Para.A2endPt-Para.A1endPt), h( (Para.A1endPt+1):(Para.A2endPt))) *... |
---|
221 | A( (1+Para.A1endPt):Para.A2endPt,:) + ... |
---|
222 | A( (1+Para.A2endPt):Para.A3endPt,:)' *... |
---|
223 | sparse(1:(Para.A3endPt-Para.A2endPt), 1:(Para.A3endPt-Para.A2endPt), h( (Para.A2endPt+1):(Para.A3endPt))) *... |
---|
224 | A( (1+Para.A2endPt):Para.A3endPt,:) + ... |
---|
225 | A( (1+Para.A3endPt):end,:)' *... |
---|
226 | sparse(1:(Para.A4endPt-Para.A3endPt), 1:(Para.A4endPt-Para.A3endPt), h( (Para.A3endPt+1):end)) *... |
---|
227 | A( (1+Para.A3endPt):end,:)) *... |
---|
228 | (2*alpha); |
---|
229 | if LogBarrierFlg |
---|
230 | hessphi_t = S'*sparse(1:length(h_t), 1:length(h_t), h_t)*S; |
---|
231 | % if VERBOSE |
---|
232 | % disp('condition number of H'); |
---|
233 | % condest((hessphi+ hessphi_t)) |
---|
234 | % end |
---|
235 | dw = -(hessphi+ hessphi_t)\(gradphi + gradphi_t); |
---|
236 | else |
---|
237 | dw = -hessphi\gradphi; |
---|
238 | end |
---|
239 | end |
---|
240 | |
---|
241 | % newton decrement=========================================== |
---|
242 | %lambdasqr = full(-(gradphi + gradphi_t)'*dw); |
---|
243 | lambdasqr = full(dw'*(hessphi+ hessphi_t)'*dw); |
---|
244 | % =========================================================== |
---|
245 | |
---|
246 | if VERBOSE |
---|
247 | disp(sprintf('%4d %15.6e %15.6e %4.5e %10.2e %10.2e %6d %10.2e %6d',... |
---|
248 | ntiter,Exact_phi, phi, alpha, max(abs( logExpTerm)), s,normg, lambdasqr/2,pflg,prelres,pitr)); |
---|
249 | end |
---|
250 | |
---|
251 | history = [history [Exact_phi + logBarrierValue; phi+logBarrierValue; ... |
---|
252 | alpha; max(abs( logExpTerm)); normg; lambdasqr/2; pcgiter]]; |
---|
253 | %------------------------------------------------------------ |
---|
254 | % STOPPING CRITERION |
---|
255 | %------------------------------------------------------------ |
---|
256 | if ( normg < ABSTOL) |
---|
257 | status = 1; |
---|
258 | disp('Minimal norm( gradient) reached.'); |
---|
259 | % disp(sprintf('%d/%d',sum(abs((A2'*h)./(2*lambda))<0.5),n)); |
---|
260 | return; |
---|
261 | end |
---|
262 | if (lambdasqr/2 <= EPSILON) |
---|
263 | status = 1; |
---|
264 | disp('Minimal Newton decrement reached'); |
---|
265 | return; |
---|
266 | end |
---|
267 | % hacking stop |
---|
268 | if ntiter ~=0 |
---|
269 | if abs(history(6,end-1) - history(6,end)) < EPSILON |
---|
270 | disp('lambdasqr not improving'); |
---|
271 | break; |
---|
272 | end |
---|
273 | end |
---|
274 | |
---|
275 | pcgiter = pcgiter+pitr; |
---|
276 | %------------------------------------------------------------ |
---|
277 | % BACKTRACKING LINE SEARCH |
---|
278 | %------------------------------------------------------------ |
---|
279 | s = 1; |
---|
280 | while max( S*(w+s*dw)+q) > 0 s = BETA*s; end % first set the new w inside |
---|
281 | for lsiter = 1:MAX_LS_ITER |
---|
282 | new_w = w+s*dw; |
---|
283 | logExpTerm = alpha*(A*new_w-b); |
---|
284 | expTerm = exp(logExpTerm); |
---|
285 | expTermNeg = exp(-logExpTerm); |
---|
286 | g = sparse( 1./(1+expTermNeg) - 1./(1+expTerm) ); |
---|
287 | newgradphi = A'*g; |
---|
288 | if (norm(newgradphi)<=(1-ALPHA_LineSearch*s)*normg) break; end |
---|
289 | s = BETA*s; |
---|
290 | end |
---|
291 | if (lsiter == MAX_LS_ITER) disp('Max lineSearch iteration reached'); break; end |
---|
292 | w = new_w; |
---|
293 | |
---|
294 | % check if new_w feasible |
---|
295 | if ~all(S*w+q<=0) |
---|
296 | disp('new_w infeasible'); |
---|
297 | break; |
---|
298 | end |
---|
299 | % tighten the sigmoid approximation |
---|
300 | alpha = min(alpha*mu, ALPHA_MAX); |
---|
301 | |
---|
302 | % tighten the log barrier approximation |
---|
303 | %t = t*mu_t; |
---|
304 | end |
---|
305 | status = -1; |
---|
306 | |
---|
307 | %------------------------------------------------------------ |
---|
308 | % CALL BACK FUNCTIONS FOR PCG |
---|
309 | %------------------------------------------------------------ |
---|
310 | %function y = AXfunc(x,A,h,d,p) |
---|
311 | % y = A'*(h.*(A*x))+d.*x; |
---|
312 | |
---|
313 | %function y = Mfunc(x,A,h,d,p) |
---|
314 | % y = x.*p; |
---|
315 | return; |
---|
316 | |
---|