Permalink
Browse files

Merge pull request #408 from vene/omp_bug

FIX for issue #403: Orthogonal Matching Pursuit bug in column swapping
  • Loading branch information...
2 parents 797041a + 1718399 commit a509d4c1c7bdb16ae8e7887e2caa047a1bd92e21 @vene vene committed Oct 21, 2011
Showing with 21 additions and 6 deletions.
  1. +6 −6 sklearn/linear_model/omp.py
  2. +15 −0 sklearn/linear_model/tests/test_omp.py
@@ -63,7 +63,7 @@ def _cholesky_omp(X, y, n_nonzero_coefs, tol=None, copy_X=True):
alpha = np.dot(X.T, y)
residual = y
n_active = 0
- idx = []
+ indices = range(X.shape[1]) # keeping track of swapping
max_features = X.shape[1] if tol is not None else n_nonzero_coefs
L = np.empty((max_features, max_features), dtype=X.dtype)
@@ -84,9 +84,9 @@ def _cholesky_omp(X, y, n_nonzero_coefs, tol=None, copy_X=True):
warn(premature)
break
L[n_active, n_active] = np.sqrt(1 - v)
- idx.append(lam)
X.T[n_active], X.T[lam] = swap(X.T[n_active], X.T[lam])
alpha[n_active], alpha[lam] = alpha[lam], alpha[n_active]
+ indices[n_active], indices[lam] = indices[lam], indices[n_active]
n_active += 1
# solves LL'x = y as a composition of two triangular systems
gamma, _ = potrs(L[:n_active, :n_active], alpha[:n_active], lower=True,
@@ -98,7 +98,7 @@ def _cholesky_omp(X, y, n_nonzero_coefs, tol=None, copy_X=True):
elif n_active == max_features:
break
- return gamma, idx
+ return gamma, indices[:n_active]
def _gram_omp(Gram, Xy, n_nonzero_coefs, tol_0=None, tol=None,
@@ -152,7 +152,7 @@ def _gram_omp(Gram, Xy, n_nonzero_coefs, tol_0=None, tol=None,
nrm2, swap = linalg.get_blas_funcs(('nrm2', 'swap'), (Gram,))
potrs, = get_lapack_funcs(('potrs',), (Gram,))
- idx = []
+ indices = range(len(Gram)) # keeping track of swapping
alpha = Xy
tol_curr = tol_0
delta = 0
@@ -176,9 +176,9 @@ def _gram_omp(Gram, Xy, n_nonzero_coefs, tol_0=None, tol=None,
warn(premature)
break
L[n_active, n_active] = np.sqrt(1 - v)
- idx.append(lam)
Gram[n_active], Gram[lam] = swap(Gram[n_active], Gram[lam])
Gram.T[n_active], Gram.T[lam] = swap(Gram.T[n_active], Gram.T[lam])
+ indices[n_active], indices[lam] = indices[lam], indices[n_active]
Xy[n_active], Xy[lam] = Xy[lam], Xy[n_active]
n_active += 1
# solves LL'x = y as a composition of two triangular systems
@@ -196,7 +196,7 @@ def _gram_omp(Gram, Xy, n_nonzero_coefs, tol_0=None, tol=None,
elif n_active == max_features:
break
- return gamma, idx
+ return gamma, indices[:n_active]
def orthogonal_mp(X, y, n_nonzero_coefs=None, tol=None, precompute_gram=False,
@@ -128,3 +128,18 @@ def test_identical_regressors():
warnings.simplefilter('always')
orthogonal_mp(newX, newy, 2)
assert len(w) == 1
+
+
+def test_swapped_regressors():
+ gamma = np.zeros(n_features)
+ # X[:, 21] should be selected first, then X[:, 0] selected second,
+ # which will take X[:, 21]'s place in case the algorithm does
+ # column swapping for optimization (which is the case at the moment)
+ gamma[21] = 1.0
+ gamma[0] = 0.5
+ new_y = np.dot(X, gamma)
+ new_Xy = np.dot(X.T, new_y)
+ gamma_hat = orthogonal_mp(X, new_y, 2)
+ gamma_hat_gram = orthogonal_mp_gram(G, new_Xy, 2)
+ assert_equal(np.flatnonzero(gamma_hat), [0, 21])
+ assert_equal(np.flatnonzero(gamma_hat_gram), [0, 21])

0 comments on commit a509d4c

Please sign in to comment.