import numpy as np
import scipy.optimize as opt
import scipy.sparse as sps
import numpy.linalg as nla
import scipy.linalg as sla
import time
def nnlsm_blockpivot(A, B, is_input_prod=False, init=None):
""" Nonnegativity-constrained least squares with block principal pivoting method and column grouping
Solves min ||AX-B||_2^2 s.t. X >= 0 element-wise.
J. Kim and H. Park, Fast nonnegative matrix factorization: An active-set-like method and comparisons,
SIAM Journal on Scientific Computing,
vol. 33, no. 6, pp. 3261-3281, 2011.
Parameters
----------
A : numpy.array, shape (m,n)
B : numpy.array or scipy.sparse matrix, shape (m,k)
Optional Parameters
-------------------
is_input_prod : True/False. - If True, the A and B arguments are interpreted as
AtA and AtB, respectively. Default is False.
init: numpy.array, shape (n,k). - If provided, init is used as an initial value for the algorithm.
Default is None.
Returns
-------
X, (success, Y, num_cholesky, num_eq, num_backup)
X : numpy.array, shape (n,k) - solution
success : True/False - True if the solution is found. False if the algorithm did not terminate
due to numerical errors.
Y : numpy.array, shape (n,k) - Y = A.T * A * X - A.T * B
num_cholesky : int - the number of Cholesky factorizations needed
num_eq : int - the number of linear systems of equations needed to be solved
num_backup: int - the number of appearances of the back-up rule. See SISC paper for details.
"""
if is_input_prod:
AtA = A
AtB = B
else:
AtA = A.T.dot(A)
if sps.issparse(B):
AtB = B.T.dot(A)
AtB = AtB.T
else:
AtB = A.T.dot(B)
(n, k) = AtB.shape
MAX_ITER = n * 5
if init is not None:
PassSet = init > 0
X, num_cholesky, num_eq = normal_eq_comb(AtA, AtB, PassSet)
Y = AtA.dot(X) - AtB
else:
X = np.zeros([n, k])
Y = -AtB
PassSet = np.zeros([n, k], dtype=bool)
num_cholesky = 0
num_eq = 0
p_bar = 3
p_vec = np.zeros([k])
p_vec[:] = p_bar
ninf_vec = np.zeros([k])
ninf_vec[:] = n + 1
not_opt_set = np.logical_and(Y < 0, ~PassSet)
infea_set = np.logical_and(X < 0, PassSet)
not_good = np.sum(not_opt_set, axis=0) + np.sum(infea_set, axis=0)
not_opt_colset = not_good > 0
not_opt_cols = not_opt_colset.nonzero()[0]
big_iter = 0
num_backup = 0
success = True
while not_opt_cols.size > 0:
big_iter += 1
if MAX_ITER > 0 and big_iter > MAX_ITER:
success = False
break
cols_set1 = np.logical_and(not_opt_colset, not_good < ninf_vec)
temp1 = np.logical_and(not_opt_colset, not_good >= ninf_vec)
temp2 = p_vec >= 1
cols_set2 = np.logical_and(temp1, temp2)
cols_set3 = np.logical_and(temp1, ~temp2)
cols1 = cols_set1.nonzero()[0]
cols2 = cols_set2.nonzero()[0]
cols3 = cols_set3.nonzero()[0]
if cols1.size > 0:
p_vec[cols1] = p_bar
ninf_vec[cols1] = not_good[cols1]
true_set = np.logical_and(not_opt_set, np.tile(cols_set1, (n, 1)))
false_set = np.logical_and(infea_set, np.tile(cols_set1, (n, 1)))
PassSet[true_set] = True
PassSet[false_set] = False
if cols2.size > 0:
p_vec[cols2] = p_vec[cols2] - 1
temp_tile = np.tile(cols_set2, (n, 1))
true_set = np.logical_and(not_opt_set, temp_tile)
false_set = np.logical_and(infea_set, temp_tile)
PassSet[true_set] = True
PassSet[false_set] = False
if cols3.size > 0:
for col in cols3:
candi_set = np.logical_or(
not_opt_set[:, col], infea_set[:, col])
to_change = np.max(candi_set.nonzero()[0])
PassSet[to_change, col] = ~PassSet[to_change, col]
num_backup += 1
(X[:, not_opt_cols], temp_cholesky, temp_eq) = normal_eq_comb(
AtA, AtB[:, not_opt_cols], PassSet[:, not_opt_cols])
num_cholesky += temp_cholesky
num_eq += temp_eq
X[abs(X) < 1e-12] = 0
Y[:, not_opt_cols] = AtA.dot(X[:, not_opt_cols]) - AtB[:, not_opt_cols]
Y[abs(Y) < 1e-12] = 0
not_opt_mask = np.tile(not_opt_colset, (n, 1))
not_opt_set = np.logical_and(
np.logical_and(not_opt_mask, Y < 0), ~PassSet)
infea_set = np.logical_and(
np.logical_and(not_opt_mask, X < 0), PassSet)
not_good = np.sum(not_opt_set, axis=0) + np.sum(infea_set, axis=0)
not_opt_colset = not_good > 0
not_opt_cols = not_opt_colset.nonzero()[0]
return X, (success, Y, num_cholesky, num_eq, num_backup)
def nnlsm_activeset(A, B, overwrite=False, is_input_prod=False, init=None):
""" Nonnegativity-constrained least squares with active-set method and column grouping
Solves min ||AX-B||_2^2 s.t. X >= 0 element-wise.
Algorithm of this routine is close to the one presented in the following paper but
is different in organising inner- and outer-loops:
M. H. Van Benthem and M. R. Keenan, J. Chemometrics 2004; 18: 441-450
Parameters
----------
A : numpy.array, shape (m,n)
B : numpy.array or scipy.sparse matrix, shape (m,k)
Optional Parameters
-------------------
is_input_prod : True/False. - If True, the A and B arguments are interpreted as
AtA and AtB, respectively. Default is False.
init: numpy.array, shape (n,k). - If provided, init is used as an initial value for the algorithm.
Default is None.
Returns
-------
X, (success, Y, num_cholesky, num_eq, num_backup)
X : numpy.array, shape (n,k) - solution
success : True/False - True if the solution is found. False if the algorithm did not terminate
due to numerical errors.
Y : numpy.array, shape (n,k) - Y = A.T * A * X - A.T * B
num_cholesky : int - the number of Cholesky factorizations needed
num_eq : int - the number of linear systems of equations needed to be solved
"""
if is_input_prod:
AtA = A
AtB = B
else:
AtA = A.T.dot(A)
if sps.issparse(B):
AtB = B.T.dot(A)
AtB = AtB.T
else:
AtB = A.T.dot(B)
(n, k) = AtB.shape
MAX_ITER = n * 5
num_cholesky = 0
num_eq = 0
not_opt_set = np.ones([k], dtype=bool)
if overwrite:
X, num_cholesky, num_eq = normal_eq_comb(AtA, AtB)
PassSet = X > 0
not_opt_set = np.any(X < 0, axis=0)
elif init is not None:
X = init
X[X < 0] = 0
PassSet = X > 0
else:
X = np.zeros([n, k])
PassSet = np.zeros([n, k], dtype=bool)
Y = np.zeros([n, k])
opt_cols = (~not_opt_set).nonzero()[0]
not_opt_cols = not_opt_set.nonzero()[0]
Y[:, opt_cols] = AtA.dot(X[:, opt_cols]) - AtB[:, opt_cols]
big_iter = 0
success = True
while not_opt_cols.size > 0:
big_iter += 1
if MAX_ITER > 0 and big_iter > MAX_ITER:
success = False
break
(Z, temp_cholesky, temp_eq) = normal_eq_comb(
AtA, AtB[:, not_opt_cols], PassSet[:, not_opt_cols])
num_cholesky += temp_cholesky
num_eq += temp_eq
Z[abs(Z) < 1e-12] = 0
infea_subset = Z < 0
temp = np.any(infea_subset, axis=0)
infea_subcols = temp.nonzero()[0]
fea_subcols = (~temp).nonzero()[0]
if infea_subcols.size > 0:
infea_cols = not_opt_cols[infea_subcols]
(ix0, ix1_subsub) = infea_subset[:, infea_subcols].nonzero()
ix1_sub = infea_subcols[ix1_subsub]
ix1 = not_opt_cols[ix1_sub]
X_infea = X[(ix0, ix1)]
alpha = np.zeros([n, len(infea_subcols)])
alpha[:] = np.inf
alpha[(ix0, ix1_subsub)] = X_infea / (X_infea - Z[(ix0, ix1_sub)])
min_ix = np.argmin(alpha, axis=0)
min_vals = alpha[(min_ix, range(0, alpha.shape[1]))]
X[:, infea_cols] = X[:, infea_cols] + \
(Z[:, infea_subcols] - X[:, infea_cols]) * min_vals
X[(min_ix, infea_cols)] = 0
PassSet[(min_ix, infea_cols)] = False
elif fea_subcols.size > 0:
fea_cols = not_opt_cols[fea_subcols]
X[:, fea_cols] = Z[:, fea_subcols]
Y[:, fea_cols] = AtA.dot(X[:, fea_cols]) - AtB[:, fea_cols]
Y[abs(Y) < 1e-12] = 0
not_opt_subset = np.logical_and(
Y[:, fea_cols] < 0, ~PassSet[:, fea_cols])
new_opt_cols = fea_cols[np.all(~not_opt_subset, axis=0)]
update_cols = fea_cols[np.any(not_opt_subset, axis=0)]
if update_cols.size > 0:
val = Y[:, update_cols] * ~PassSet[:, update_cols]
min_ix = np.argmin(val, axis=0)
PassSet[(min_ix, update_cols)] = True
not_opt_set[new_opt_cols] = False
not_opt_cols = not_opt_set.nonzero()[0]
return X, (success, Y, num_cholesky, num_eq)
def normal_eq_comb(AtA, AtB, PassSet=None):
""" Solve many systems of linear equations using combinatorial grouping.
M. H. Van Benthem and M. R. Keenan, J. Chemometrics 2004; 18: 441-450
Parameters
----------
AtA : numpy.array, shape (n,n)
AtB : numpy.array, shape (n,k)
Returns
-------
(Z,num_cholesky,num_eq)
Z : numpy.array, shape (n,k) - solution
num_cholesky : int - the number of unique cholesky decompositions done
num_eq: int - the number of systems of linear equations solved
"""
num_cholesky = 0
num_eq = 0
if AtB.size == 0:
Z = np.zeros([])
elif (PassSet is None) or np.all(PassSet):
Z = nla.solve(AtA, AtB)
num_cholesky = 1
num_eq = AtB.shape[1]
else:
Z = np.zeros(AtB.shape)
if PassSet.shape[1] == 1:
if np.any(PassSet):
cols = PassSet.nonzero()[0]
Z[cols] = nla.solve(AtA[np.ix_(cols, cols)], AtB[cols])
num_cholesky = 1
num_eq = 1
else:
#
# Both _column_group_loop() and _column_group_recursive() work well.
# Based on preliminary testing,
# _column_group_loop() is slightly faster for tiny k(<10), but
# _column_group_recursive() is faster for large k's.
#
grps = _column_group_recursive(PassSet)
for gr in grps:
cols = PassSet[:, gr[0]].nonzero()[0]
if cols.size > 0:
ix1 = np.ix_(cols, gr)
ix2 = np.ix_(cols, cols)
#
# scipy.linalg.cho_solve can be used instead of numpy.linalg.solve.
# For small n(<200), numpy.linalg.solve appears faster, whereas
# for large n(>500), scipy.linalg.cho_solve appears faster.
# Usage example of scipy.linalg.cho_solve:
# Z[ix1] = sla.cho_solve(sla.cho_factor(AtA[ix2]),AtB[ix1])
#
Z[ix1] = nla.solve(AtA[ix2], AtB[ix1])
num_cholesky += 1
num_eq += len(gr)
num_eq += len(gr)
return Z, num_cholesky, num_eq
def _column_group_loop(B):
""" Given a binary matrix, find groups of the same columns
with a looping strategy
Parameters
----------
B : numpy.array, True/False in each element
Returns
-------
A list of arrays - each array contain indices of columns that are the same.
"""
initial = [np.arange(0, B.shape[1])]
before = initial
after = []
for i in range(0, B.shape[0]):
all_ones = True
vec = B[i]
for cols in before:
if len(cols) == 1:
after.append(cols)
else:
all_ones = False
subvec = vec[cols]
trues = subvec.nonzero()[0]
falses = (~subvec).nonzero()[0]
if trues.size > 0:
after.append(cols[trues])
if falses.size > 0:
after.append(cols[falses])
before = after
after = []
if all_ones:
break
return before
def _column_group_recursive(B):
""" Given a binary matrix, find groups of the same columns
with a recursive strategy
Parameters
----------
B : numpy.array, True/False in each element
Returns
-------
A list of arrays - each array contain indices of columns that are the same.
"""
initial = np.arange(0, B.shape[1])
return [a for a in column_group_sub(B, 0, initial) if len(a) > 0]
def column_group_sub(B, i, cols):
vec = B[i][cols]
if len(cols) <= 1:
return [cols]
if i == (B.shape[0] - 1):
col_trues = cols[vec.nonzero()[0]]
col_falses = cols[(~vec).nonzero()[0]]
return [col_trues, col_falses]
else:
col_trues = cols[vec.nonzero()[0]]
col_falses = cols[(~vec).nonzero()[0]]
after = column_group_sub(B, i + 1, col_trues)
after.extend(column_group_sub(B, i + 1, col_falses))
return after
def _test_column_grouping(m=10, n=5000, num_repeat=5, verbose=False):
print ('\nTesting column_grouping ...\n')
A = np.array([[True, False, False, False, False],
[True, True, False, True, True]])
grps1 = _column_group_loop(A)
grps2 = _column_group_recursive(A)
grps3 = [np.array([0]),
np.array([1, 3, 4]),
np.array([2])]
print ('OK' if all([np.array_equal(a, b) for (a, b) in zip(grps1, grps2)]) else 'Fail')
print ('OK' if all([np.array_equal(a, b) for (a, b) in zip(grps1, grps3)]) else 'Fail')
for i in iter(range(0, num_repeat)):
A = np.random.rand(m, n)
B = A > 0.5
start = time.time()
grps1 = _column_group_loop(B)
elapsed_loop = time.time() - start
start = time.time()
grps2 = _column_group_recursive(B)
elapsed_recursive = time.time() - start
if verbose:
print ('Loop :', elapsed_loop)
print ('Recursive:', elapsed_recursive)
print ('OK' if all([np.array_equal(a, b) for (a, b) in zip(grps1, grps2)]) else 'Fail')
# sorted_idx = np.concatenate(grps)
# print B
# print sorted_idx
# print B[:,sorted_idx]
return
def _test_normal_eq_comb(m=10, k=3, num_repeat=5):
print ('\nTesting normal_eq_comb() ...\n')
for i in iter(range(0, num_repeat)):
A = np.random.rand(2 * m, m)
X = np.random.rand(m, k)
C = (np.random.rand(m, k) > 0.5)
X[~C] = 0
B = A.dot(X)
B = A.T.dot(B)
A = A.T.dot(A)
Sol, a, b = normal_eq_comb(A, B, C)
print ('OK' if np.allclose(X, Sol) else 'Fail')
return
def _test_nnlsm():
print ('\nTesting nnls routines ...\n')
m = 100
n = 10
k = 200
rep = 5
for r in iter(range(0, rep)):
A = np.random.rand(m, n)
X_org = np.random.rand(n, k)
X_org[np.random.rand(n, k) < 0.5] = 0
B = A.dot(X_org)
# B = np.random.rand(m,k)
# A = np.random.rand(m,n/2)
# A = np.concatenate((A,A),axis=1)
# A = A + np.random.rand(m,n)*0.01
# B = np.random.rand(m,k)
import time
start = time.time()
C1, info = nnlsm_blockpivot(A, B)
elapsed2 = time.time() - start
rel_norm2 = nla.norm(C1 - X_org) / nla.norm(X_org)
print ('nnlsm_blockpivot: ', 'OK ' if info[0] else 'Fail',\
'elapsed:{0:.4f} error:{1:.4e}'.format(elapsed2, rel_norm2))
start = time.time()
C2, info = nnlsm_activeset(A, B)
num_backup = 0
elapsed1 = time.time() - start
rel_norm1 = nla.norm(C2 - X_org) / nla.norm(X_org)
print ('nnlsm_activeset: ', 'OK ' if info[0] else 'Fail',\
'elapsed:{0:.4f} error:{1:.4e}'.format(elapsed1, rel_norm1))
import scipy.optimize as opt
start = time.time()
C3 = np.zeros([n, k])
for i in iter(range(0, k)):
res = opt.nnls(A, B[:, i])
C3[:, i] = res[0]
elapsed3 = time.time() - start
rel_norm3 = nla.norm(C3 - X_org) / nla.norm(X_org)
print ('scipy.optimize.nnls: ', 'OK ',\
'elapsed:{0:.4f} error:{1:.4e}'.format(elapsed3, rel_norm3))
if num_backup > 0:
break
if rel_norm1 > 10e-5 or rel_norm2 > 10e-5 or rel_norm3 > 10e-5:
break
print ('')
if __name__ == '__main__':
_test_column_grouping()
_test_normal_eq_comb()
_test_nnlsm()
