[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, alfa, status, history, T_nt_hist] = ... |
---|
| 40 | SigmoidSolver( Para, t_0, alpha_0, w_0, 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) |
---|
| 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 | % |
---|
| 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 | |
---|
| 87 | global A b |
---|
| 88 | |
---|
| 89 | |
---|
| 90 | % FEASIBLAE DATA CHECK |
---|
| 91 | [m,n] = size(A); % problem size: m examples, n features |
---|
| 92 | |
---|
| 93 | %------------------------------------------------------------ |
---|
| 94 | % INITIALIZE |
---|
| 95 | %------------------------------------------------------------ |
---|
| 96 | |
---|
| 97 | % SIGMOID APPROXIMATION |
---|
| 98 | NUMERICAL_LIMIT_EXP = 3e2; % this value is calculated such that exp(2*NUMERICAL_LIMIT_EXP) < inf |
---|
| 99 | MU_alpha = 1; ;%1.1; %1.5; % for sigmoidal |
---|
| 100 | MU_alpha2 = 1.3; %5;%1.5; %1.5; % for sigmoidal |
---|
| 101 | ALPHA_MAX = 50000; |
---|
| 102 | 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 |
---|
| 103 | ALPHA_MIN = 5000.0; % should be atleast 100 -- ERROR |
---|
| 104 | |
---|
| 105 | % NEWTON PARAMETERS |
---|
| 106 | MAX_TNT_ITER = 1; %100; %25; % maximum Newton iteration |
---|
| 107 | %ABSTOL = 1e-8; % terminates when the norm of gradient < ABSTOL |
---|
| 108 | EPSILON = 1e-6; %1e-6; % terminate when lambdasqr_by_2 < EPSILON |
---|
| 109 | |
---|
| 110 | % LINE SEARCH PARAMETERS |
---|
| 111 | ALPHA_LineSearch= 0.01; % minimum fraction of decrease in norm(gradient) |
---|
| 112 | BETA = 0.5; % stepsize decrease factor |
---|
| 113 | MAX_LS_ITER = 25; % maximum backtracking line search iteration |
---|
| 114 | %MIN_STEP_SIZE = 1e-20; |
---|
| 115 | |
---|
| 116 | |
---|
| 117 | % VARIABLE INITIALZE |
---|
| 118 | w = zeros(size(A,2),1); |
---|
| 119 | |
---|
| 120 | %save initial_w.mat w |
---|
| 121 | % load initial_w.mat |
---|
| 122 | |
---|
| 123 | % setting starting alpha_0 |
---|
| 124 | alpha_0 = 1 / max( abs(A*w-b) ); |
---|
| 125 | alfa = alpha_0; |
---|
| 126 | alfa_max = alfa; |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | % ----------------------------------------------------------- |
---|
| 130 | % META LOOP FOR CHANGING t and ALPHA |
---|
| 131 | % ---------------------------------------------------------- |
---|
| 132 | history=[]; |
---|
| 133 | T_nt_hist=[]; |
---|
| 134 | |
---|
| 135 | if VERBOSE |
---|
| 136 | disp(sprintf('%15s %15s %11s %12s %10s %10s %7s %15s %15s',... |
---|
| 137 | 'exact obj', 'primal obj', 'alpha', 'MlogExpTerm',... |
---|
| 138 | 'norm(g)', 'lambda^2By2', ... |
---|
| 139 | 'gap', 'search_step', 'newton_steps')); |
---|
| 140 | %'normg_sigmoid', 'normg_t', ... |
---|
| 141 | end |
---|
| 142 | |
---|
| 143 | status = -1; |
---|
| 144 | %for outer_iter = 1:MAX_OUTER_ITER |
---|
| 145 | %history_Inner = []; |
---|
| 146 | % Aw_b = ; |
---|
| 147 | logExpTerm = alfa./MU_alpha2*(A*w-b)'; |
---|
| 148 | expTerm = exp(logExpTerm); |
---|
| 149 | dw = zeros(n,1); % dw newton step |
---|
| 150 | while alfa < ALPHA_MIN |
---|
| 151 | |
---|
| 152 | % history_Inner = []; |
---|
| 153 | |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | % Makes the profiler not work! |
---|
| 157 | % if (VERBOSE >= 2) |
---|
| 158 | % disp(sprintf('%s %15s %15s %11s %12s %10s %10s %10s %10s %9s %6s %6s %10s %10s %10s',... |
---|
| 159 | % 'iter', 'exact obj', 'primal obj', 'alpha', 'MlogExpTerm',... |
---|
| 160 | % 'In_step', 'stepsize','norm(g)', 'norm(dw)', ... |
---|
| 161 | % 'lambda^2By2','normg_sigmoid','normg_t', ... |
---|
| 162 | % 'normdw_t', 'normdw_sigmoid')); |
---|
| 163 | % end |
---|
| 164 | |
---|
| 165 | %------------------------------------------------------------ |
---|
| 166 | % MAIN LOOP |
---|
| 167 | %------------------------------------------------------------ |
---|
| 168 | %total_linesearch = 0; |
---|
| 169 | %ntiter = 0; |
---|
| 170 | %initialize Newton step |
---|
| 171 | |
---|
| 172 | logExpTerm = logExpTerm.*MU_alpha2; |
---|
| 173 | expTerm = expTerm.^MU_alpha2; |
---|
| 174 | %expTermNeg = exp(-logExpTerm); |
---|
| 175 | expTerm_Rec = 1./(1+expTerm); |
---|
| 176 | |
---|
| 177 | expTermNeg_Rec = 1./(1+exp(-logExpTerm)); |
---|
| 178 | g_sigmoid = (expTermNeg_Rec - expTerm_Rec); |
---|
| 179 | gradphi = (g_sigmoid*A)'; %A1endPt) + ... |
---|
| 180 | % gradphi = gradphi_sigmoid; |
---|
| 181 | |
---|
| 182 | % newtonLoop = (ntiter <= MAX_TNT_ITER); |
---|
| 183 | % while newtonLoop |
---|
| 184 | %ntiter = ntiter + 1; |
---|
| 185 | |
---|
| 186 | indicesInf = find( (logExpTerm > NUMERICAL_LIMIT_EXP) | (logExpTerm < -NUMERICAL_LIMIT_EXP) ); |
---|
| 187 | %indicesOkay = setdiff(1:length(logExpTerm), [indicesInf_Pos, indicesInf_Neg]); |
---|
| 188 | |
---|
| 189 | h_sigmoid = expTerm .* (expTerm_Rec.^2); |
---|
| 190 | h_sigmoid(indicesInf) = 0; |
---|
| 191 | |
---|
| 192 | |
---|
| 193 | %------------------------------------------------------------ |
---|
| 194 | % CALCULATE NEWTON STEP |
---|
| 195 | %------------------------------------------------------------ |
---|
| 196 | |
---|
| 197 | hessphi = (sparse(1:m,1:m,h_sigmoid)*A)' * A; |
---|
| 198 | % hessphi_sigmoid = (h_sigmoid*(A.^2))'; %diagonal element only |
---|
| 199 | |
---|
| 200 | % if DEBUG_FLAG && (any(h_sigmoid<0) || any(h_t<0) ) |
---|
| 201 | % disp(sprintf('h negative = %5.4e; h_t negative = %5.4e',h_sigmoid(h_sigmoid<0), h_t(h_t<0))); |
---|
| 202 | % end |
---|
| 203 | |
---|
| 204 | |
---|
| 205 | % hessphi_t = (h_t*(S.^2))' / (2*alfa); %Diagonal Element only |
---|
| 206 | |
---|
| 207 | %hessphi = hessphi_sigmoid; %HACK |
---|
| 208 | |
---|
| 209 | dw = - hessphi \ gradphi; |
---|
| 210 | |
---|
| 211 | % dw = - (1./hessphi) .* gradphi; % Diagonal Hessian |
---|
| 212 | dw = dw / (2*alfa); |
---|
| 213 | |
---|
| 214 | % newton decrement=========================================== |
---|
| 215 | % lambdasqr = full(-( gradphi'*dw) ); |
---|
| 216 | %lambdasqr2 = full(dw'*(hessphi+ hessphi_t)'*dw); |
---|
| 217 | % =========================================================== |
---|
| 218 | |
---|
| 219 | |
---|
| 220 | %------------------------------------------------------------ |
---|
| 221 | % BACKTRACKING LINE SEARCH |
---|
| 222 | %------------------------------------------------------------ |
---|
| 223 | s = 1; |
---|
| 224 | normg = norm(gradphi); |
---|
| 225 | lsiter = 0; |
---|
| 226 | backIterationLoop = true; |
---|
| 227 | while backIterationLoop |
---|
| 228 | lsiter = lsiter+1; |
---|
| 229 | new_w = w+s*dw; |
---|
| 230 | % Aw_b = (A*new_w-b); |
---|
| 231 | logExpTerm = alfa*(A*new_w-b)'; |
---|
| 232 | expTerm = exp(logExpTerm); |
---|
| 233 | |
---|
| 234 | expTerm_Rec = 1./(1+expTerm); |
---|
| 235 | expTermNeg_Rec = 1./(1+exp(-logExpTerm)); |
---|
| 236 | g_sigmoid = (expTermNeg_Rec - expTerm_Rec); |
---|
| 237 | gradphi = (g_sigmoid*A)'; %A1endPt) + ... |
---|
| 238 | %gradphi = gradphi_sigmoid; |
---|
| 239 | |
---|
| 240 | backIterationLoop = (lsiter <= MAX_LS_ITER) && ( norm(gradphi) > (1-ALPHA_LineSearch*s)*normg); |
---|
| 241 | s = BETA*s; |
---|
| 242 | end |
---|
| 243 | |
---|
| 244 | % total_linesearch = total_linesearch + lsiter; |
---|
| 245 | |
---|
| 246 | % if VERBOSE >= 2 |
---|
| 247 | % Exact_phi_sigmoid = norm( logExpTerm/alfa, 1); % only for debug |
---|
| 248 | % Exact_f = Exact_phi_sigmoid + logBarrierValue / t; % MIN -- fixed |
---|
| 249 | % %normdw = norm(dw); |
---|
| 250 | % normg_sigmoid = norm(gradphi_sigmoid); |
---|
| 251 | % normg_t = norm(gradphi_t); |
---|
| 252 | % normdw_sigmoid = 0;%norm( dw_sigmoid); |
---|
| 253 | % normdw_t = 0;%norm( dw_t); |
---|
| 254 | % normdw = norm( dw); |
---|
| 255 | % f_new = phi_sigmoid + logBarrierValue / t; % MIN -- fixed ? |
---|
| 256 | % disp(sprintf('%4d %15.6e %15.6e %4.5e %10.2e %10.2e %10.2e %6d %6d %10.2e %6d %6d %6d %6d',... |
---|
| 257 | % ntiter,Exact_f, f_new, alfa, max(abs( logExpTerm)), s0, s, normg, normdw, lambdasqr/2, ... |
---|
| 258 | % normg_sigmoid, normg_t, normdw_t, normdw_sigmoid)); |
---|
| 259 | % end |
---|
| 260 | |
---|
| 261 | %------------------------------------------------------------ |
---|
| 262 | % STOPPING CRITERION |
---|
| 263 | %------------------------------------------------------------ |
---|
| 264 | |
---|
| 265 | % if (lambdasqr/2 <= EPSILON) |
---|
| 266 | % status = 1; |
---|
| 267 | % %disp('Minimal Newton decrement reached'); |
---|
| 268 | % end |
---|
| 269 | |
---|
| 270 | % newtonLoop = (ntiter <= MAX_TNT_ITER) && (lsiter < MAX_LS_ITER); % |
---|
| 271 | |
---|
| 272 | % set new w |
---|
| 273 | w = new_w; |
---|
| 274 | %alfa = min(alfa*MU_alpha, alfa_max); |
---|
| 275 | % If you want to do the above, i.e., increase the value in each |
---|
| 276 | % Newton step, uncomment the gradient calculation, if ntiter == 1 |
---|
| 277 | |
---|
| 278 | % end % -----------END of the MAIN LOOP----------------------------- |
---|
| 279 | |
---|
| 280 | % Tighten the sigmoid and log approximation |
---|
| 281 | |
---|
| 282 | |
---|
| 283 | % if VERBOSE |
---|
| 284 | % Exact_phi_sigmoid = norm( logExpTerm/alfa, 1); % only for debug |
---|
| 285 | % |
---|
| 286 | % normg_sigmoid = norm(gradphi_sigmoid); |
---|
| 287 | % normdw = norm( dw); |
---|
| 288 | % normg = norm( gradphi); |
---|
| 289 | % expTermNeg = exp(-logExpTerm); |
---|
| 290 | % phi_sigmoid = ( log( 1+ expTermNeg ) + log( 1+ expTerm) ); |
---|
| 291 | % phi_sigmoid( indicesInf_Pos ) = logExpTerm( indicesInf_Pos ); |
---|
| 292 | % phi_sigmoid( indicesInf_Neg ) = -logExpTerm( indicesInf_Neg ); |
---|
| 293 | % phi_sigmoid = (1/alfa) * sum(phi_sigmoid,2); |
---|
| 294 | % |
---|
| 295 | % f_new = phi_sigmoid; % MIN -- fixed ? |
---|
| 296 | % Exact_f = Exact_phi_sigmoid; % MIN -- fixed |
---|
| 297 | % disp(sprintf('%15.6e %15.6e %4.5e %10.2e %10.2e %10.2e %10.2d %10d %10d',... |
---|
| 298 | % Exact_f, f_new, alfa, max(abs( logExpTerm)), normg, lambdasqr/2, ... |
---|
| 299 | % gap, total_linesearch, ntiter)); |
---|
| 300 | % % normg_sigmoid, normg_t, ... |
---|
| 301 | % end |
---|
| 302 | |
---|
| 303 | % T_nt_hist = [T_nt_hist history_Inner]; |
---|
| 304 | % history=[history [length(history_Inner)]]; |
---|
| 305 | |
---|
| 306 | %if alfa >= ALPHA_MIN/MU_alpha |
---|
| 307 | %else |
---|
| 308 | % epsilon_gap = (1-alfa/ALPHA_MIN)*1e2 + (alfa/ALPHA_MIN)*EPSILON_GAP; |
---|
| 309 | %end |
---|
| 310 | |
---|
| 311 | |
---|
| 312 | % alfa = min(alfa*MU_alpha, ALPHA_MAX); |
---|
| 313 | % alfa_max = min(alfa_max*MU_alpha2, ALPHA_MAX); |
---|
| 314 | % if (alfa > ALPHA_MIN) && (status >= 1) |
---|
| 315 | % t = MU_t*t; % if just waiting for gap, then increase it at double speed |
---|
| 316 | % end |
---|
| 317 | |
---|
| 318 | alfa = alfa*MU_alpha2; |
---|
| 319 | % alfa_max = min(alfa*MU_alpha2, ALPHA_MAX); |
---|
| 320 | |
---|
| 321 | end |
---|
| 322 | |
---|
| 323 | return; |
---|
| 324 | |
---|