Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

WIP: Covariance Updates #911

Open
wants to merge 37 commits into from

3 participants

ibayer Alexandre Gramfort Olivier Grisel
ibayer

This very much work in process for issue #910.

The idea is to get first a python version of the covariant updates to work and turn it then step by
step into cython code.

Reference:
[1] Friedman, J., T. Hastie, and R. Tibshirani. “Regularization Paths for Generalized Linear Models via Coordinate Descent.” Journal of Statistical Software 33, no. 1 (2010): 1.

ibayer

I though I could hijack tmp from the current implementation as this seems to be the z variable of the soft-thresholding operator in Eq. 6 . The resulting code is now to slow to test on real data but since test_enet_toy() is failing something must be wrong. I don't know if this is due to an coding error or a conceptual problem. I can't find the 1/N term from Eq. 8 in the current implementation?

Alexandre Gramfort
Owner

isn't why your trying to do what enet_coordinate_descent_gram is doing ie work with a precomputed gram matrix XTX and XTy ?

ibayer

Not exactly I want to precompute only XTy. Components of XTX are only computed on demand "each time a new feature X[:,k], enters the model (for the first time), we need to compute and store its inner product with all the rest of the features" [1]. But it's pretty similar to enet_coordinate_descent_gram and I wasn't aware of that. Thanks

ibayer

I'm stuck can't get naive implementation of Eq. 9 working in enet_coordinate_descent . :( I will try to get there by modifying enet_coordinate_descent_gram next. equations I assume that z in equation (6) is tmp in the current implementation. Any suggestions are welcome.

Alexandre Gramfort
Owner

yes z is tmp in cd_fast.pyx

try writing your own pure python code so you can debug. You can quickly reimplement enet_coordinate_descent in pure python to see that you understand it (forget the dual gap for now) and start from there.

ibayer

I have implemented the glmnet algorithm in python unit-testing every line. I get the same coef as with the modified cd_fast.pxy for X = [-1, 0, 1], y = [-1, 0, 1]. Doing all calculation by hand I get the same coef too.

I did a visual local convergence check, checking where the calculated coef are on the enet objective function. I found that the the difference between the calculated coef and the optimal is noticeable bigger then the difference in the objective function. The result is in the neighborhood of the optimal point but not on it.

I did some more experiments with X element of R^n. I found that the results are very close to the optimum after the soft-thresholding operator but the following "proportional shrinkage for the rige penalty" (divisor in Eq (5)) moves the coef away from the optimum. Something else that I realized is that the algorithm converges after one step, but I think that should be expected.

I can't explain the observed behavior. I must have an conceptual error somewhere do you have an idea? I will next reimplement enet_coordinate_descent line by line as suggested maybe that helps me to find the missing piece.

Alexandre Gramfort
Owner
ibayer

Reimplementing and right debugging did the trick, thanks fore the valuable advice!
I think the following tricked my for quite some time (in coordinate_descent.py):

alpha = self.alpha * self.rho * n_samples
beta = self.alpha * (1.0 - self.rho) * n_samples

I would suggest renaming them (at least alpha) it's kind a mean to explain what alpha is and then change
it in subsequent code. Here the meaning of parameters alpha and beta are not clear as well sparse_enet_coordinate_descent(np.ndarray[DOUBLE, ndim=1] w,
double alpha, double beta ... )
` .

I'll put up code soon.

Alexandre Gramfort
Owner
ibayer

One Test was not working properly, fixed it now it dosn't pass anymore :( .
https://gist.github.com/2978601
For the Naive Updates test_2d passes only when the accuracy is limited to 4 digits and test_big_data 2 digits. The objective function decreases for 5 steps and then jumps up and down keeping the first 6 digits constant. Is this a bug or numerical noise?
The Covariance Updates are commented out, they fail completely.

Alexandre Gramfort
Owner
ibayer

Thanks a lot, this was starting to drive me crazy :) .
I just had to change tol to get the tests passing. Could you please look at the Covariance Updates, I checked the calculations by hand but just can't find the problem.

It should be a straight forward implementation of Eq. 8 replacing the fist term on the right hand side by Eq. 9.

I suspect the 1/N to cause the problem. The 1/N devision is done implicitly and has to be redone for the beta_j term in Eq 8.

sklearn/linear_model/cd_fast.pyx
@@ -148,6 +148,12 @@ def enet_coordinate_descent(np.ndarray[DOUBLE, ndim=1] w,
<DOUBLE*>(X.data + ii * n_samples * sizeof(DOUBLE)), 1,
<DOUBLE*>R.data, 1)
+ # test code for Covariance Updates
+ # right hand side of Eq. 8, the sum is replaced by the right hand side
+ # of Eq. 9
+ tmp_sum = np.array([np.dot(X[:, ii], X[:, k]) * w[k] for k in xrange(n_features)]).sum()
+ tmp = (1/n_features)*(np.dot(X[:, ii], y) - tmp_sum) + w_ii
Alexandre Gramfort Owner

1/n_features will be equal to 0 if n_features is an int which it should be

Alexandre Gramfort Owner

use 1. / n_features or 1 / float(n_features) to avoid this

Alexandre Gramfort Owner

tell me if it fixes the problem

ibayer
ibayer added a note

I'm sorry I should have been more explicit about which code I'm talking.
The int/int wasn't the problem in this case, but thanks for the reminder I still fall into
this trap from time to time.
I was referring to the out commented code in your gist:
https://gist.github.com/2979853

Alexandre Gramfort Owner
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
ibayer

Covariance Updates work now but are still without gradient caching and updating.
https://gist.github.com/2984543

Alexandre Gramfort
Owner
ibayer

Feature inner product caching and gradient updates are working.
Iterate only over active features, is next.
https://gist.github.com/2984543

Alexandre Gramfort
Owner
ibayer

Iteration only over active features works in separate implementation:
https://gist.github.com/2984543
But test_lasso_path() fails if the same is included in cd_fast.pyx with "Objective did not converge", all other tests pass. The same behavior can be observed if one includes:

            if n_iter > 0 and w[ii] == 0:
                continue

to skip updates for features that are not active after the first pass.

Alexandre Gramfort
Owner
ibayer

I already increased from 150 to 3000. But the test should not start failing if I skip iterations that are not expected to change ẁ anyway. All other tests kept passing ( test_enet_path() etc. )! I'll keep debugging but I'm somewhat in the dark.

Alexandre Gramfort
Owner
ibayer

I'm now using a user specified memory limit for the caching of the feature_inner_product . I don't use dynamic memory allocation but check after each outer iteration if the number of active features is small enough for caching. If that's the case caching is used from that point on. I want to do some refactoring and then start profiling. Especially the active set storing and updating need some thinking.

Alexandre Gramfort
Owner

I am following.

let me cc @paolo-losi so he can also give us a hand reviewing.

are you aware of "cython -a" to identify bottlenecks in cython code?

Olivier Grisel
Owner

What profiler are use using? Do you know yep and google-profiler? http://pypi.python.org/pypi/yep

There is also valgrind and http://www.vrplumber.com/programming/runsnakerun/ with meliae for memory analysis.

ibayer

@ogrisel I tried to use yep and google-profiler but the output makes not much sense to me.
Output and profile script:
https://gist.github.com/3027998
I get the following output if I run the profile script:
PROFILE: interrupts/evictions/bytes = 0/0/32
I don't know if this is expected or what is means.
I didn't find much information on how to use google-profiler in general.

The implementation of the active set of features could be improved.
Even though the algorithm only loop the active features (coordinates) all gradients
get updated

gradient -= feature_inner_product[ii, :] * (w[ii] - w_ii)

but if would be sufficient to do something like:

gradient[active_set] -= feature_inner_product[ii, active_set] *   (w[ii] - w_ii)

but the indexing costs more then the saved operations. This is also true for other matrix/vector operations.
@agramfort I'm going to look into "cython -a" now but this only gives information about potential bottlenecks (as I understand it).

Alexandre Gramfort
Owner
ibayer

I haven't found a clever way around it yet. Shrinking the vectors before calling np.dot is expensive as well, calling cblas directly makes the indexing complicated. Using a loop over the active indexes will probably loose even against larger blas operations.

ibayer

I just realized that "You removed the blas calls so you have no chance for now to reach the speed of the previous implementation." was a general comment. Isn't np.dot basically a call to blas dot? If it speeds things up I will switch to daxpy etc. but I want to do this as late as possible as it complicates debugging.

Alexandre Gramfort
Owner
Alexandre Gramfort
Owner
ibayer

Okay I understand. Does the out factoring in cython functions, as I did with the dual gap calculation, impact performance too? I was hopping to make things a bit easier to read.

Alexandre Gramfort
Owner
Alexandre Gramfort
Owner

rather than fixing a size limit for gram you can investigate some dynamic memory allocation strategy using np.resize

it's not a big deal if you use np.resize in cython as long as you don't do it too often.

ibayer

Thanks I didn't know about np.resize I will look into it. I added the last test to test a bug I just found, I have to fix this first. I was thinking about some semi dynamic memory allocation strategy too, like re-size every time half the entries in w are zero. But I'm not sure how much speedup this really brings. If the main time is spend without caching it won't matter much. If the function is called in a path algorithm fashion a resizing will happen anyway every time the function is called. I need to type def everything and use blas to be able to do performance evaluations.

Alexandre Gramfort
Owner
Olivier Grisel
Owner

I think a set of cython utilities to handle resizable numpy arrays would be very interesting and a candidate for reuse in other parts of the project.

sklearn/linear_model/cd_fast.pyx
@@ -83,12 +83,58 @@ def sparse_std(unsigned int n_samples,
@cython.boundscheck(False)
@cython.wraparound(False)
+cdef inline double calculate_gap(np.ndarray[DOUBLE, ndim=1] w,
+ double alpha, double beta,
+ np.ndarray[DOUBLE, ndim=2] X,
+ np.ndarray[DOUBLE, ndim=1] y, bint positive):
ibayer
ibayer added a note

I keep getting

warning: cd_fast.pyx:86:33: Buffer unpacking not optimized away.
warning: cd_fast.pyx:86:33: Buffer unpacking not optimized away.
warning: cd_fast.pyx:88:28: Buffer unpacking not optimized away.
warning: cd_fast.pyx:88:28: Buffer unpacking not optimized away.
warning: cd_fast.pyx:89:28: Buffer unpacking not optimized away.
warning: cd_fast.pyx:89:28: Buffer unpacking not optimized away.

Something seems to be wrong with my function declarations, could someone give me a hint?

Olivier Grisel Owner
ogrisel added a note

I think we could start using the memory view API of cython 0.16.

http://docs.cython.org/src/userguide/memoryviews.html

Unless of course you need to call numpy specific methods that have no memory views syntax equivalent.

Alternatively you can pass array pointers using the & operator as I did here:

dbea9d8

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
ibayer

The one python call left in the inner loop, is for array resizing:

v = np.array([0, 4, 6., 0, 0, 22, 1, 99., 0, 0])
index = np.nonzero(v)[0]
v = v[index] # = np.resize(v, len(index))

I'm trying to do this in cython but without success yet. Any suggestions?

Olivier Grisel
Owner

Allocate a new array and copy the content of the previous array into it. If possible it would be better to implement an exponentially auto-growing array based datastructure as @agramfort told previously.

ibayer

I want to shrink the array, I was hoping this would be possible without reallocation. Copying is straight forward but I don't know how to do the shrinking without relaying on python functions as show above. I already tried something like this:

def shrink_vector(np.ndarray[DOUBLE, ndim=1] v,
                  np.ndarray[INTEGER, ndim=1] index, int length):
    cdef int i
    for i in range(length):
        v[i] = v[index[i]]

but that involves python operations as well. (That's at least how I understand the cython -a output)

Olivier Grisel
Owner

Actually, w.r.t. my previous comments, there is already an ArrayBuilder utility in scikit-learn (implemented by @larsmans):

https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/utils/arraybuilder.pyx

It used in the cython code for the sklearn.datasets.svmlight_format package. However it does not support shrinking.

For you use case, the shuffling implemented inplace is probably false

cdef int i
    for i in range(length):
        v[i] = v[index[i]]

For instance if index = [1, 0] and v = [42, 43] you will get v = [43, 43] which is probably not what you are expecting. You probably need a (possibly preallocated) tmp array to perform the shuffling in a safe manner.

Olivier Grisel ogrisel closed this
Olivier Grisel ogrisel reopened this
Olivier Grisel
Owner

Sorry wrong button.

ibayer

@ogrisel thanks for the tip, but ArrayBuilder uses np.resize() too.
The following should work:

def shrink_vector(np.ndarray[DOUBLE, ndim=1] v,
                  np.ndarray[INTEGER, ndim=1] index, int length):
    cdef int i
    for i in range(length):
        v[i] = v[index[i]]
    v = np.resize(v, length)
    return v

I can do to the copying in place but the np.resize() makes a new copy and is a python function. I was hoping to avoid both. I will try to get the cython profiling working, maybe all this doesn't matter much anyway.

Alexandre Gramfort
Owner
Olivier Grisel
Owner

I am also wondering if the following signature should be a cdef to avoid being seen as a python function callable from the module:

+1. If need it both internal and as a public utility you can also use the cpdef keyword instead.

ibayer

I need to profile the cython code, but all that yep gave me is:

Using local file /home/mane/virt_env/scikit-learn/bin/python.
Using local file speed_test.py.prof.
Total: 15 samples
       2  13.3%  13.3%        2  13.3% __open
       1   6.7%  80.0%        1   6.7% PyArray_DescrFromType
       0   0.0%  80.0%        1   6.7% PyObject_Malloc
       1   6.7%  86.7%        1   6.7% PyObject_Malloc (inline)
       1   6.7%  93.3%        1   6.7% __xstat64
       1   6.7% 100.0%        1   6.7% _dl_rtld_di_serinfo

if I call this script:

import cd_fast as fast
import numpy as np
from sklearn.datasets.samples_generator import make_regression

X, y = make_regression(n_samples=200, n_features=5, n_informative=5,
                        random_state=0)
rho = 0.80
alpha = 5

n_samples, n_features = X.shape
l2_reg = alpha * rho * n_samples
l1_reg = alpha * (1.0 - rho) * n_samples
w = np.zeros(n_features)
X = np.asfortranarray(X)
result, gap, tol = fast.enet_coordinate_descent(w, l2_reg, l1_reg,
                                X, y, max_iter=100, tol=1e-7, positive=False)

with python -m yep -v speed_test.py. I spend quit some time looking for other cython profiling solutions...

Olivier Grisel
Owner

15 samples is very short to draw any meaningful conclusion from this profiling session. How long does it take to run your script ? You should probably grow the dataset enough so that the script runs in ~30s.

Olivier Grisel
Owner

BTW, there is no such think as a cython profiler: python runs a compiled extension that is written in C (that was generated from the cython source). In the output of the profiler you will merely see symbols from the compiled extension which is to be expected but should be enough to understand where the bottleneck is.

ibayer

@ogrisel 15 samples is from the yep output, I assume it means that the script has been executed 15 times.
""In the output of the profiler you will merely see symbols from the compiled extension which is to be expected but should be enough to understand where the bottleneck is.""
Makes sense, but a cython profiler would have been nice :( . Thanks for mention it, I'll look for information on how to profile c extensions now.

Olivier Grisel
Owner

No yep does launch the script only once IIRC. I don't know the default sampling frequency of the google perftools sampling profiler. How long does an individual run last ?

You should make it last enough for the sampling profiler to sample enough stacktraces to compute interesting statistics on the interesting region (the part that does the glm optimization).

ibayer

< 5 sec , I'm going to increase it.

%timeit result, gap, tol = fast.enet_coordinate_descent(w, l2_reg, l1_reg,X, y, max_iter=100, tol=1e-7, positive=False, memory_limit=1000000)
1 loops, best of 3: 7.77 s per loop
Olivier Grisel
Owner

I did some tests and apparently the dataset generation part is dominating. Try to add those lines near the beginning.

from sklearn.externals import joblib
memory = joblib.Memory('.')
make_regression = memory.cache(make_regression)

Then run it once (e.g. with n_samples=1000 and n_features=10000) for building the cache and the second time for profiling.

Olivier Grisel
Owner

Also to read the output:

Analyzing Text Output

Text mode has lines of output that look like this:

       14   2.1%  17.2%       58   8.7% std::_Rb_tree::find

Here is how to interpret the columns:

    Number of profiling samples in this function
    Percentage of profiling samples in this function
    Percentage of profiling samples in the functions printed so far
    Number of profiling samples in this function and its callees
    Percentage of profiling samples in this function and its callees
    Function name 

from: http://google-perftools.googlecode.com/svn/trunk/doc/cpuprofile.html

ibayer

"Then run it once (e.g. with n_samples=1000 and n_features=10000) for building the cache and the second time for profiling."
Running via "python -m yep -v speed_test.py"? (can't see how the cache could be preserved)
Thanks

Olivier Grisel
Owner

joblib will use the local folder to store the generated dataset.

sklearn/linear_model/cd_fast.pyx
((27 lines not shown))
+ if (dual_norm_XtA > alpha):
+ const = alpha / dual_norm_XtA
+ A_norm = R_norm * const
+ gap = 0.5 * (R_norm ** 2 + A_norm ** 2)
+ else:
+ const = 1.0
+ gap = R_norm ** 2
+
+ gap += alpha * linalg.norm(w, 1) - const * np.dot(R.T, y) + \
+ 0.5 * beta * (1 + const ** 2) * (w_norm ** 2)
+ return gap
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef inline restore_w(np.ndarray[DOUBLE, ndim=1] w,
Alexandre Gramfort Owner

you can use cpdef on this function as it's private

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
ibayer

"You can essentially think of a cpdef method as a cdef method + some extras." http://docs.cython.org/src/userguide/pyrex_differences.html I don't seed the advantage to make it cpdef, do you think we need to call it from Python?
Some profiling results:
with caching:
result, gap, tol = fast.enet_coordinate_descent(w, l2_reg, l1_reg,
X, y, max_iter=10000, tol=1e-7, positive=False, memory_limit=100000000)

Total: 3888 samples
    2665  68.5%  68.5%     2665  68.5% DOUBLE_dot
     272   7.0%  75.5%      272   7.0% ATL_dgemvT_a1_x1_b0_y1
     180   4.6%  80.2%      180   4.6% ATL_daxpy_xp1yp1aXbX
      77   2.0%  82.2%       77   2.0% PyArray_MatrixProduct2
      54   1.4%  83.5%       54   1.4% _aligned_contig_to_strided_size8
       1   0.0%  83.6%       49   1.3% iter_ass_subscript
      48   1.2%  84.8%       48   1.2% iter_ass_sub_int (inline)
      36   0.9%  85.7%       45   1.2% DOUBLE_nonzero
      39   1.0%  86.7%       39   1.0% DOUBLE_absolute
      33   0.8%  87.6%       35   0.9% DOUBLE_copyswap
      32   0.8%  88.4%       32   0.8% getaliasbyname_r
      23   0.6%  89.0%       27   0.7% __pyx_pf_7sklearn_12linear_model_7cd_fast_2enet_coordinate_descent
      26   0.7%  89.7%       26   0.7% _IO_str_pbackfail
      25   0.6%  90.3%       25   0.6% ATL_ddot_xp0yp0aXbX
      17   0.4%  90.7%       17   0.4% PyArray_Nonzero

The first line DOUBLE_dot really suprises me, AFAIK does this mean that most of the time is spend in a numpy.dot function. Dataset generation: X, y = make_regression(n_samples=100, n_features=5000, n_informative=30, random_state=0) . I'm just reading how I can get more information from google-profiler.

Alexandre Gramfort
Owner
Olivier Grisel
Owner

Use the joblib cache as I said previously to skip the dataset generation in the profiling script (run it once without profiling to fill in the cache and the second time with yep to actually run the profiler).

Olivier Grisel
Owner

Alternatively use the yep.start() / yep.stop() calls around the elastic net optimization line in to restrict profiling to that portion of the script: http://pypi.python.org/pypi/yep

ibayer

@ogrisel that's how I generated the profiling output above. Using the yep.start() / yep.stop() gives the same result if I issue top20 during an interactive google-pprof session. Is there something in the output that made you believe otherwise?
@agramfort I was using a high alpha but increasing it again, it made only a minor difference.

rho = 0.80
alpha = 50
ibayer

I think the high number of samples in DOUBLE_dot comes from checking the duality gap. After out-commenting the dual gap call. This is the new output:

Total: 565 samples
     296  52.4%  52.4%      296  52.4% ATL_dgemvT_a1_x1_b0_y1
      52   9.2%  61.6%       52   9.2% _aligned_contig_to_strided_size8
      47   8.3%  69.9%       50   8.8% DOUBLE_nonzero
      34   6.0%  75.9%       34   6.0% ATL_ddot_xp0yp0aXbX
      15   2.7%  78.6%       15   2.7% iter_subscript_int (inline)
      14   2.5%  81.1%       14   2.5% __memset_sse2_rep
      13   2.3%  83.4%       16   2.8% __pyx_pf_7sklearn_12linear_model_7cd_fast

The current implementation does not show this behaviour, I'm therefore investigating why the dual gap check is called so often.

Olivier Grisel
Owner

This is indeed what I was suspecting. I think we should have an option to disable the duality gap convergence check when doing regularization path optimiziation (and just rely upon the normalized absolute difference on w as a convergence check). I think this is glmnet does and seems much less expensive (this needs to be checked).

ibayer

What I don't understand is why the dual gap checking is slowing down my glmnet implementation but not the current Naive Updates. Is this expected behaviour or is it a bug? The glmnet convergence check is even simpler:

"If another complete cycle does not change the active set, we are done, otherwise the process is repeated. Active-set convergence is also mentioned in Meier et al. [2008] and Krishnapuram and Hartemink [2005]" [1]

sklearn/linear_model/cd_fast.pyx
((8 lines not shown))
- if alpha == 0:
- warnings.warn("Coordinate descent with alpha=0 may lead to unexpected"
+ if l2_reg == 0:
+ warnings.warn("Coordinate descent with l2_reg=0 may lead to unexpected"
" results and is discouraged.")
Alexandre Gramfort Owner

should be l1_reg=0

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
sklearn/linear_model/cd_fast.pyx
((103 lines not shown))
if positive and tmp < 0:
w[ii] = 0.0
else:
- w[ii] = fsign(tmp) * fmax(fabs(tmp) - alpha, 0) \
- / (norm_cols_X[ii] + beta)
+ w[ii] = fsign(tmp) * fmax(fabs(tmp) - l2_reg, 0) \
Alexandre Gramfort Owner

you inverted l1_reg and l2_reg

alpha is the coefficient that penalizes the L1 norm

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Alexandre Gramfort
Owner
ibayer

@agramfort "Then do a full pass on all features." This would be very expensive, especially if the active set is small and the cache of not active features is released. Why not using the dual gap check only if the active set doesn't change anymore?
I would agree with @ogrisel that dual gap checks or any other expensive convergence test makes little sense if enet_coordinate_descent is used as part of a regularization path optimization. It would be sufficient to check the dual gap at the last point on the path.

Alexandre Gramfort
Owner

@agramfort "Then do a full pass on all features." This would be very expensive, especially if the active set is small

that's what the current code is doing once for each iteration and it's not particularly slow.

and the cache of not active features is released.

the full pass can be done with no covariance precomputation.

Why not using the dual gap check only if the active set doesn't change anymore?

why not but if dual_gap is not used as a convergence check than the semantic of tol
will change (just to keep in mind).

Olivier Grisel
Owner

The meaning of tol should be documented. Also it might be interesting to make it possible to enable the dual convergence check for each path step as a boolean option to make it possible for the user to trade speed for accuracy.

ibayer

While hunting down the reason why the dual gap check is called so often in the active set implementation I found the following:
(without removal of zeros coef)
In a particular case a coefficient, that has been zero for some iterations, ends up with a value of -0.000942165848007.
The dual gap for this run is up to 6 digits equal to the dual gap of the current implementation.
(with removal)
If I use the active set strategy this coefficient is removed from the active set and ends up with a value of 0. The dual gap is here 45.4270949122. All coefficients appear to be equal up to 4 digits. I used tol=1e-7.

I assume that the high dual gap value is responsible for the many dual gap checks in this case.

Alexandre Gramfort
Owner
ibayer
Alexandre Gramfort
Owner
ibayer
Alexandre Gramfort
Owner

tell me when I can review. Glad you're moving forward !

Alexandre Gramfort
Owner

here is gist that @fabianp wrote this morning to help you out:

https://gist.github.com/3097107

it's not the covariance update but the implementation of strong rules.

ibayer

@fabianp
thanks a lot, I have a closer look as soon I'm done with re-factoring my covariance updates with active-set iterations.
@agramfort
I think a came up with a efficient solution for the active-set iterations. I'm not shrinking and growing w and the other vectors as the active set changes but order the vectors so that the active features are left aligned in the array. This has the advantage that I can simply use the first x position to do the cblas operations only on the active-set. Growing and shrinking of the vectors can then be realized by simple swaps. This works already but all the indexing makes the code hard to follow, I'm currently re-factoring everything. I'll update you asas things are in a stable state again.

Alexandre Gramfort
Owner
ibayer

Hey, can't come up with a good answer to the following:
I'm swapping random elements in a vector, how can I restore the original order?
What information do I need to store about the swaps?

Alexandre Gramfort
Owner
ibayer

Thanks, I think that's similar to what I was doing.
But if I'm randomly swapping elements and one element can be involved
in more then one swap, then it's not clear to me how to build this permutation
vector. The permutation vector from step i get's permuted in step i+1. So the
updating of the permutation vector is the main problem.

I was recording the swaps by doing the same swaps to index = np.arange(len(a)) but I think that breaks down as soon as one element is involved more then once. I'll think it over.

Alexandre Gramfort
Owner
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Jun 17, 2012
  1. ibayer

    slow python prototype

    ibayer authored
Commits on Jun 18, 2012
  1. ibayer

    still not working

    ibayer authored
Commits on Jun 26, 2012
  1. ibayer

    Covariance Updates added

    ibayer authored
  2. ibayer

    old code removed

    ibayer authored
Commits on Jun 27, 2012
  1. ibayer

    active set iteration added

    ibayer authored
Commits on Jun 28, 2012
  1. ibayer

    add old enet_coordinate_descent for regression tests,

    ibayer authored
    factor out dual gap calculation
  2. ibayer

    factor out active set reduction

    ibayer authored
  3. ibayer

    user specified memory limit added. If the limit is to low, values are

    ibayer authored
    recomputed and user warning issued.
  4. ibayer
Commits on Jun 29, 2012
  1. ibayer
  2. ibayer

    testing cd_fast directly

    ibayer authored
Commits on Jun 30, 2012
  1. ibayer

    gradient updates vectorized

    ibayer authored
  2. ibayer

    avoid loop, by not making a difference between initialized and not

    ibayer authored
    initialized gradients. The values will be overwritten during
    initialization anyway.
  3. ibayer
  4. ibayer
Commits on Jul 1, 2012
  1. ibayer

    add some type defs

    ibayer authored
Commits on Jul 2, 2012
  1. ibayer

    Add test to check if cache indexing works if cache is not sufficiently

    ibayer authored
    large to use it for all features.
Commits on Jul 3, 2012
  1. ibayer

    use <= instead of <

    ibayer authored
  2. ibayer

    fix indexing error

    ibayer authored
  3. ibayer
Commits on Jul 4, 2012
  1. ibayer

    some cblas operations added

    ibayer authored
Commits on Jul 5, 2012
  1. ibayer

    cblas operations added

    ibayer authored
  2. ibayer

    cblas dgemv added

    ibayer authored
  3. ibayer

    feature_inner_product change from fortran to C continuous, this makes

    ibayer authored
    cached iterations faster but not cached iterations slower.
  4. ibayer
Commits on Jul 9, 2012
  1. ibayer
Commits on Jul 10, 2012
  1. ibayer
Commits on Jul 11, 2012
  1. ibayer
  2. ibayer
  3. ibayer
  4. ibayer
  5. ibayer
  6. ibayer

    simplifications

    ibayer authored
Commits on Jul 12, 2012
  1. ibayer

    simplifications

    ibayer authored
  2. ibayer

    make indexing more explicit

    ibayer authored
  3. ibayer

    simplifications

    ibayer authored
  4. ibayer

    fix wrong index

    ibayer authored
Something went wrong with that request. Please try again.