MNT: Use nrm2 to validate and normalize atoms in dictionary learning
@@ -8,7 +8,7 @@ import time
import sys
import itertools
-from math import sqrt, ceil
+from math import ceil
import numpy as np
from scipy import linalg
@@ -378,6 +378,7 @@ def _update_dict(dictionary, Y, code, verbose=False, return_r2=False,
# Get BLAS functions
gemm, = linalg.get_blas_funcs(('gemm',), (dictionary, code, Y))
ger, = linalg.get_blas_funcs(('ger',), (dictionary, code))
+ nrm2, = linalg.get_blas_funcs(('nrm2',), (dictionary,))
# Residuals, computed with BLAS for speed and efficiency
# R <- -1.0 * U * V^T + 1.0 * Y
# Outputs R as Fortran array for efficiency
@@ -389,8 +390,9 @@ def _update_dict(dictionary, Y, code, verbose=False, return_r2=False,
if positive:
np.clip(dictionary[:, k], 0, None, out=dictionary[:, k])
# Scale k'th atom
- atom_norm_square = np.dot(dictionary[:, k], dictionary[:, k])
- if atom_norm_square < 1e-20:
+ # (U_k * U_k) ** 0.5
+ atom_norm = nrm2(dictionary[:, k])
+ if atom_norm < 1e-10:
if verbose == 1:
sys.stdout.write("+")
sys.stdout.flush()
@@ -401,10 +403,11 @@ def _update_dict(dictionary, Y, code, verbose=False, return_r2=False,
np.clip(dictionary[:, k], 0, None, out=dictionary[:, k])
# Setting corresponding coefs to 0
code[k, :] = 0.0
- dictionary[:, k] /= sqrt(np.dot(dictionary[:, k],
- dictionary[:, k]))
+ # (U_k * U_k) ** 0.5
+ atom_norm = nrm2(dictionary[:, k])
+ dictionary[:, k] /= atom_norm
else:
- dictionary[:, k] /= sqrt(atom_norm_square)
+ dictionary[:, k] /= atom_norm
# R <- -1.0 * U_k * V_k^T + R
R = ger(-1.0, dictionary[:, k], code[k, :], a=R, overwrite_a=True)
if return_r2:
