[37] | 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 | |
---|