Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 59 additions & 20 deletions lib/scikit/decomposition.flow
Original file line number Diff line number Diff line change
Expand Up @@ -765,6 +765,32 @@ export function fast_ica_free(model: FastICA) -> void {
# SparsePCA
# ============================================================================

# Convergence statistic shared by sparse_pca_fit and minibatch_sparse_pca_fit
# (issue #460). Returns the largest absolute change in any coordinate of the
# component across one full update pass: gradient step, soft-threshold,
# renormalize.
#
# Max rather than mean: the solution of an L1-penalized objective is sparse, so
# most coordinates sit at exactly zero and stop moving early. A mean over p
# would divide that handful of still-moving active coordinates by the full
# feature count and report convergence while the support is still shifting. The
# max is what the iteration actually has to drive to zero.
#
# Both arguments are live buffers owned by the caller; nothing here reads freed
# memory (issue #378).
#
# Module-unique name per compiler bug #465.
function _spca_max_abs_delta(prev: ptr<f32>, cur: ptr<f32>, p: i32) -> f32 {
let mut worst: f32 = 0.0
for j in 0 to p {
let d: f32 = fabs((cur[j] - prev[j]) as f64) as f32
if d > worst {
worst = d
}
}
return worst
}

export struct SparsePCA {
components: ptr<ptr<f32> >,
n_components: i32,
Expand Down Expand Up @@ -816,9 +842,15 @@ export function sparse_pca_fit(X: Matrix, n_components: i32, alpha: f32, max_ite
}
}

let prev: ptr<f32> = array_new_f32(p)

for iter in 0 to max_iter {
n_iter = iter + 1

for j in 0 to p {
prev[j] = components[c][j]
}

let grad: ptr<f32> = array_new_f32(p)
for i in 0 to n {
let mut dot: f32 = 0.0
Expand Down Expand Up @@ -860,19 +892,20 @@ export function sparse_pca_fit(X: Matrix, n_components: i32, alpha: f32, max_ite
}
}

if iter > 0 && iter % 10 == 0 {
let mut change: f32 = 0.0
for j in 0 to p {
change = change + fabs((grad[j]) as f64) as f32
}
if change / (p as f32) < tol {
array_free_f32(grad)
break
}
}
# Issue #460: stop on the change in the parameters, not on the
# gradient magnitude. The gradient of an L1 objective does not
# vanish at a soft-thresholded solution, so the old test could
# never reach tol and every fit ran the full max_iter.
let delta: f32 = _spca_max_abs_delta(prev, components[c], p)

array_free_f32(grad)

if delta < tol {
break
}
}

array_free_f32(prev)
}

matrix_free(Xc)
Expand Down Expand Up @@ -969,9 +1002,15 @@ export function minibatch_sparse_pca_fit(
}
}

let prev: ptr<f32> = array_new_f32(p)

for iter in 0 to max_iter {
n_iter = iter + 1

for j in 0 to p {
prev[j] = components[c][j]
}

# Sample mini-batch
let actual_batch: i32 = batch_size
if actual_batch > n { actual_batch = n }
Expand Down Expand Up @@ -1020,19 +1059,19 @@ export function minibatch_sparse_pca_fit(
}
}

if iter > 0 && iter % 10 == 0 {
let mut change: f32 = 0.0
for j in 0 to p {
change = change + fabs((grad[j]) as f64) as f32
}
if change / (p as f32) < tol {
array_free_f32(grad)
break
}
}
# Issue #460: same parameter-change criterion as sparse_pca_fit.
# The mini-batch gradient is noisier still, so a gradient-magnitude
# test is even less meaningful here.
let delta: f32 = _spca_max_abs_delta(prev, components[c], p)

array_free_f32(grad)

if delta < tol {
break
}
}

array_free_f32(prev)
}

matrix_free(Xc)
Expand Down
31 changes: 24 additions & 7 deletions tests/test_new_modules_v2.flow
Original file line number Diff line number Diff line change
Expand Up @@ -886,19 +886,36 @@ function test_sparse_pca() -> i32 {
println(n_zeros)
return 1
}
# The convergence test compares the raw gradient magnitude against tol,
# which never falls below 0.0001 on this workload, so the loop runs the
# full 200 iterations. Under the bug it stopped at iteration 10 whenever
# the freed buffer happened to read back as near-zero.
if spca.n_iter != 200 {
print(" FAIL: expected n_iter 200, got ")
# Issue #460 replaced the convergence test. It used to compare the raw
# gradient magnitude against tol, which an L1 objective never drives to
# zero, so the loop always ran the full 200 iterations and tol did
# nothing. It now compares the largest absolute change in the component
# over one update pass, which the soft-thresholded iteration does drive to
# zero. The loadings are identical to the ones the 200-iteration run
# produced; only the iteration count moved.
#
# Two accepted values, because n_iter here is platform-dependent. The
# component initialisation consumes libc rand(), and Darwin's Park-Miller
# generator gives a different starting vector from glibc's TYPE_3
# generator, so the fit reaches its fixed point at a different iteration.
# A standalone C replica of this fixture and of sparse_pca_fit, checked to
# reproduce the Flow output exactly on Darwin, gives 19 under Darwin libc
# and 24 under glibc on both linux/amd64 and linux/arm64. The zero count
# is 6 on every one of them.
#
# This assertion is still the issue #378 gate: it pins an exact iteration
# count, so a convergence test reading freed memory would make it vary
# from run to run on one unmodified binary.
if spca.n_iter != 19 && spca.n_iter != 24 {
print(" FAIL: expected n_iter 19 (Darwin) or 24 (glibc), got ")
println(spca.n_iter)
return 1
}

print(" OK: 2 components, ")
print(n_zeros)
println(" sparse zeros in loading vectors")
print(" sparse zeros in loading vectors, n_iter = ")
println(spca.n_iter)

matrix_free(transformed)
sparse_pca_free(spca)
Expand Down
Loading
Loading