Skip to content

Commit cbc5418

Browse files
committed
WIP
1 parent 7c9fbe1 commit cbc5418

File tree

2 files changed

+42
-2
lines changed

2 files changed

+42
-2
lines changed

skglm/datafits/single_task.py

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -158,11 +158,10 @@ def initialize(self, X, y):
158158
def initialize_sparse(self, X_data, X_indptr, X_indices, y):
159159
n_features = len(X_indptr) - 1
160160
self.lipschitz = np.zeros(n_features, dtype=X_data.dtype)
161-
self.global_lipschitz = 0.
162161
for j in range(n_features):
163162
Xj = X_data[X_indptr[j]:X_indptr[j+1]]
164163
self.lipschitz[j] = (Xj ** 2).sum() / (len(y) * 4)
165-
self.global_lipschitz += (Xj ** 2).sum() / (len(y) * 4)
164+
self.global_lipschitz = 0.
166165

167166
def value(self, y, w, Xw):
168167
return np.log(1. + np.exp(- y * Xw)).sum() / len(y)

skglm/utils.py

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -486,3 +486,44 @@ def prox_vec(w, z, penalty, lipschitz):
486486
for j in range(n_features):
487487
w[j] = penalty.prox_1d(z[j], 1 / lipschitz, j)
488488
return w
489+
490+
491+
@njit
492+
def power_method(X_data, X_indptr, X_indices, n_iter):
493+
"""Power method to compute largest eigenvalue of sparse matrix X.
494+
495+
Parameters
496+
----------
497+
X_data : array, shape (n_elements,)
498+
`data` attribute of the sparse CSC matrix X.
499+
500+
X_indptr : array, shape (n_features + 1,)
501+
`indptr` attribute of the sparse CSC matrix X.
502+
503+
X_indices : array, shape (n_elements,)
504+
`indices` attribute of the sparse CSC matrix X.
505+
506+
n_iter : int
507+
Number of iterations for the power method.
508+
509+
Returns
510+
-------
511+
rayleigh : float
512+
Rayleigh quotient or spectral radius of X^TX
513+
"""
514+
n_features = len(X_indptr) - 1
515+
b_k = np.random.rand(n_features)
516+
for _ in range(n_iter):
517+
b_k1 = np.zeros(n_features)
518+
for j in range(n_features):
519+
b_k1[j] =
520+
b_k1 = A @ b_k # write the full loop
521+
b_k1_nrm = norm(b_k1)
522+
b_k = b_k1 / b_k1_nrm
523+
rayleigh = b_k1 @ b_k / (norm(b_k) ** 2)
524+
return rayleigh
525+
526+
527+
528+
529+

0 commit comments

Comments
 (0)