#include "mex.h" #include using namespace std; // Dynamic programming algorithm for the 1d fused lasso problem // (Ryan's implementation of Nick Johnson's algorithm) // ported to matlab // % before calling this function, please run // % \$ mex prox_matlab.cpp // % in matlab to compile. void prox_dp(int n, double *y, double lam, double *beta) { // Take care of a few trivial cases if (n==0) return; if (n==1 || lam==0) { for (int i=0; i -lam) break; alo += a[lo]; blo += b[lo]; } // Compute the negative knot tm[k] = (-lam-blo)/alo; l = lo-1; x[l] = tm[k]; // Compute hi: step down from r until the // derivative is less than lam ahi = alast; bhi = blast; for (hi=r; hi>=l; hi--) { if (-ahi*x[hi]-bhi < lam) break; ahi += a[hi]; bhi += b[hi]; } // Compute the positive knot tp[k] = (lam+bhi)/(-ahi); r = hi+1; x[r] = tp[k]; // Update a and b a[l] = alo; b[l] = blo+lam; a[r] = ahi; b[r] = bhi+lam; afirst = 1; bfirst = -lam-y[k+1]; alast = -1; blast = -lam+y[k+1]; } // Compute the last coefficient: this is where // the function has zero derivative alo = afirst; blo = bfirst; for (lo=l; lo<=r; lo++) { if (alo*x[lo]+blo > 0) break; alo += a[lo]; blo += b[lo]; } beta[n-1] = -blo/alo; // Compute the rest of the coefficients, by the // back-pointers for (int k=n-2; k>=0; k--) { if (beta[k+1]>tp[k]) beta[k] = tp[k]; else if (beta[k+1]