From 5f48e74c3864ff960130e86d7be7b75b611dba79 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 17 May 2020 11:35:13 +0200 Subject: [PATCH] Enforced dtype of input matrix and data is mantained Changes in spgl1 and some auxiliary routines are required to ensure that the dtype of the solution is consistent with that of the input matrix and data. This is mostly done by specifying the correct dtype when creating new arrays. --- pytests/test_spgl1.py | 137 +++++++++++++++++++++++++++++++----------- spgl1/spgl1.py | 14 ++--- 2 files changed, 110 insertions(+), 41 deletions(-) diff --git a/pytests/test_spgl1.py b/pytests/test_spgl1.py index ca2c823..64b8bfb 100755 --- a/pytests/test_spgl1.py +++ b/pytests/test_spgl1.py @@ -3,7 +3,7 @@ from numpy.testing import assert_array_almost_equal from scipy.sparse import spdiags, csr_matrix -from spgl1.spgl1 import norm_l12nn_primal, norm_l12nn_dual, \ +from spgl1.spgl1 import oneprojector, norm_l12nn_primal, norm_l12nn_dual, \ norm_l12nn_project from spgl1 import spg_lasso, spg_bp, spg_bpdn, spg_mmv @@ -11,33 +11,57 @@ from spgl1.lsqr import lsqr -# dense matrix +# dense matrix (dtype=np.float64) par1 = {'n': 50, 'm': 50, 'k': 14, - 'sparse': False} # square + 'dtype': np.float64, 'sparse': False} # square par2 = {'n': 128, 'm': 50, 'k': 14, - 'sparse': False} # underdetermined + 'dtype': np.float64, 'sparse': False} # underdetermined par3 = {'n': 50, 'm': 100, 'k': 14, - 'sparse': False} # overdetermined + 'dtype': np.float64, 'sparse': False} # overdetermined -# sparse matrix +# sparse matrix (dtype=np.float64) par1_sp = {'n': 50, 'm': 50, 'k': 14, - 'sparse': True} # square + 'dtype': np.float64, 'sparse': True} # square par2_sp = {'n': 128, 'm': 50, 'k': 14, - 'sparse': True} # underdetermined + 'dtype': np.float64, 'sparse': True} # underdetermined par3_sp = {'n': 50, 'm': 100, 'k': 14, - 'sparse': True} # overdetermined + 'dtype': np.float64, 'sparse': True} # overdetermined -np.random.seed(1) +# dense matrix (dtype=np.float32) +par1_f32 = {'n': 50, 'm': 50, 'k': 14, + 'dtype': np.float32, 'sparse': False} # square +par2_f32 = {'n': 128, 'm': 50, 'k': 14, + 'dtype': np.float32, 'sparse': False} # underdetermined +par3_f32 = {'n': 50, 'm': 100, 'k': 14, + 'dtype': np.float32, 'sparse': False} # overdetermined + + +@pytest.mark.parametrize("par", [(par1), (par1_f32)]) +def test_oneprojector(par): + """One projector with tau=1 + + minimize_x ||b-x||_2 subject to ||x||_1 <= 1. + """ + np.random.seed(0) + + b = np.random.normal(0., 1., par['n']).astype(dtype=par['dtype']) + d = 1. + bp = oneprojector(b, d, 1.) + assert bp.dtype == par['dtype'] + assert np.linalg.norm(bp, 1) < 1. @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_lasso(par): """LASSO problem for ||x||_1 <= pi: minimize ||Ax-b||_2 subject to ||x||_1 <= 3.14159... """ + np.random.seed(0) + n = par['n'] m = par['m'] k = par['k'] @@ -47,11 +71,12 @@ def test_lasso(par): A = A1.copy() if par['sparse']: A = csr_matrix(A) + A = A.astype(par['dtype']) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m) + x = np.zeros(m, dtype=par['dtype']) x[p] = np.random.randn(k) # Set up vector b, and run solver @@ -59,21 +84,26 @@ def test_lasso(par): tau = np.pi xinv, resid, _, _ = spg_lasso(A, b, tau, verbosity=0) - assert np.linalg.norm(xinv, 1) - np.pi < 1e-10 + assert xinv.dtype == par['dtype'] + assert np.linalg.norm(xinv, 1) - np.pi < 1e-5 # Run solver with subspace minimization xinv, resid, _, _ = spg_lasso(A, b, tau, subspace_min=True, verbosity=0) - assert np.linalg.norm(xinv, 1) - np.pi < 1e-10 + assert xinv.dtype == par['dtype'] + assert np.linalg.norm(xinv, 1) - np.pi < 1e-5 @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_bp(par): """Basis pursuit (BP) problem: minimize ||x||_1 subject to Ax = b """ + np.random.seed(0) + # Create random m-by-n encoding matrix n = par['n'] m = par['m'] @@ -82,31 +112,37 @@ def test_bp(par): A, A1 = np.linalg.qr(np.random.randn(n, m), 'reduced') if m > n: A = A1.copy() + A = A.astype(par['dtype']) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m) + x = np.zeros(m, dtype=par['dtype']) x[p] = np.random.randn(k) # Set up vector b, and run solver b = A.dot(x) xinv, _, _, _ = spg_bp(A, b, verbosity=0) + assert xinv.dtype == par['dtype'] assert_array_almost_equal(x, xinv, decimal=3) # Run solver with subspace minimization xinv, _, _, _ = spg_bp(A, b, subspace_min=True, verbosity=0) + assert xinv.dtype == par['dtype'] assert_array_almost_equal(x, xinv, decimal=3) @pytest.mark.parametrize("par", [(par1), (par3), - (par1_sp), (par3_sp)]) + (par1_sp), (par3_sp), + (par1_f32), (par3_f32)]) def test_bpdn(par): """Basis pursuit denoise (BPDN) problem: minimize ||x||_1 subject to ||Ax - b||_2 <= 0.1 """ + np.random.seed(0) + # Create random m-by-n encoding matrix n = par['n'] m = par['m'] @@ -115,32 +151,42 @@ def test_bpdn(par): A, A1 = np.linalg.qr(np.random.randn(n, m), 'reduced') if m > n: A = A1.copy() + A = A.astype(par['dtype']) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m) + x = np.zeros(m, dtype=par['dtype']) x[p] = np.random.randn(k) # Set up vector b, and run solver - b = A.dot(x) + np.random.randn(n) * 0.075 + b = A.dot(x) + np.random.randn(n).astype(par['dtype']) * 0.075 + sigma = 0.10 xinv, resid, _, _ = spg_bpdn(A, b, sigma, iter_lim=1000, verbosity=0) - assert np.linalg.norm(resid) < sigma*1.1 # need to check why resid is slighly bigger than sigma + assert xinv.dtype == par['dtype'] + assert np.linalg.norm(resid) < sigma * 1.1 # need to check why resid is slighly bigger than sigma # Run solver with subspace minimization xinv, _, _, _ = spg_bpdn(A, b, sigma, subspace_min=True, verbosity=0) - assert np.linalg.norm(resid) < sigma*1.1 # need to check why resid is slighly bigger than sigma + assert xinv.dtype == par['dtype'] + assert np.linalg.norm(resid) < sigma * 1.1 # need to check why resid is slighly bigger than sigma @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_bp_complex(par): """Basis pursuit (BP) problem for complex variables: minimize ||x||_1 subject to Ax = b """ + np.random.seed(0) + + dtypecomplex = (np.ones(1, dtype=par['dtype']) + + 1j * np.ones(1, dtype=par['dtype'])).dtype + # Create random m-by-n encoding matrix n = par['n'] m = par['m'] @@ -149,25 +195,29 @@ def test_bp_complex(par): A, A1 = np.linalg.qr(np.random.randn(n, m), 'reduced') if m > n: A = A1.copy() + A = A.astype(dtypecomplex) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m, dtype=np.complex) + x = np.zeros(m, dtype=dtypecomplex) x[p] = np.random.randn(k) + 1j * np.random.randn(k) # Set up vector b, and run solver b = A.dot(x) xinv, _, _, _ = spg_bp(A, b, verbosity=0) + assert xinv.dtype == dtypecomplex assert_array_almost_equal(x, xinv, decimal=3) # Run solver with subspace minimization xinv, _, _, _ = spg_bp(A, b, subspace_min=True, verbosity=0) + assert xinv.dtype == dtypecomplex assert_array_almost_equal(x, xinv, decimal=3) @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_weighted_bp(par): """Weighted Basis pursuit (WBP) problem: @@ -178,6 +228,8 @@ def test_weighted_bp(par): minimize ||Wx||_1 subject to Ax = b """ + np.random.seed(0) + # Create random m-by-n encoding matrix n = par['n'] m = par['m'] @@ -186,31 +238,36 @@ def test_weighted_bp(par): A, A1 = np.linalg.qr(np.random.randn(n, m), 'reduced') if m > n: A = A1.copy() + A = A.astype(par['dtype']) # Create sparse vector p = np.random.permutation(m) p = p[0:k] - x = np.zeros(m) + x = np.zeros(m, dtype=par['dtype']) x[p] = np.random.randn(k) # Set up weights w and vector b - w = 0.1*np.random.rand(m) + 0.1 # Weights + w = 0.1 * np.random.rand(m).astype(dtype=par['dtype']) + 0.1 # Weights b = A.dot(x / w) # Signal xinv, _, _, _ = spg_bp(A, b, iter_lim=1000, weights=w, verbosity=0) # Reconstructed solution, with weighting xinv *= w + assert xinv.dtype == par['dtype'] assert_array_almost_equal(x, xinv, decimal=3) @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_multiple_measurements(par): """Multiple measurement vector (MMV) problem: minimize ||Y||_1,2 subject to AW^{-1}Y = B """ + np.random.seed(0) + # Create random m-by-n encoding matrix m = par['m'] n = par['n'] @@ -218,12 +275,14 @@ def test_multiple_measurements(par): l = 6 A = np.random.randn(m, n) + A = A.astype(par['dtype']) + p = np.random.permutation(n)[:k] - X = np.zeros((n, l)) + X = np.zeros((n, l), dtype=par['dtype']) X[p, :] = np.random.randn(k, l) - weights = 0.1 * np.random.rand(n) + 0.1 - W = 1 / weights * np.eye(n) + weights = 0.1 * np.random.rand(n).astype(par['dtype']) + 0.1 + W = 1 / weights * np.eye(n, dtype=par['dtype']) B = A.dot(W).dot(X) @@ -234,15 +293,20 @@ def test_multiple_measurements(par): X_w, _, _, _ = spg_mmv(A, B, 0, weights=weights, verbosity=0) X_w = spdiags(weights, 0, n, n).dot(X_w) - assert_array_almost_equal(X, X_uw, decimal=2) - assert_array_almost_equal(X, X_w, decimal=2) + assert X_uw.dtype == par['dtype'] + assert X_w.dtype == par['dtype'] + assert_array_almost_equal(X, X_uw, decimal=1) + assert_array_almost_equal(X, X_w, decimal=1) @pytest.mark.parametrize("par", [(par1), (par2), (par3), - (par1_sp), (par2_sp), (par3_sp)]) + (par1_sp), (par2_sp), (par3_sp), + (par1_f32), (par2_f32), (par3_f32)]) def test_multiple_measurements_nonnegative(par): """Multiple measurement vector (MMV) problem with non-negative norm: """ + np.random.seed(0) + # Create random m-by-n encoding matrix m = par['m'] n = par['n'] @@ -250,8 +314,10 @@ def test_multiple_measurements_nonnegative(par): l = 6 A = np.random.randn(m, n) + A = A.astype(par['dtype']) + p = np.random.permutation(n)[:k] - X = np.zeros((n, l)) + X = np.zeros((n, l), dtype=par['dtype']) X[p, :] = np.abs(np.random.randn(k, l)) B = A.dot(X) @@ -260,6 +326,7 @@ def test_multiple_measurements_nonnegative(par): primal_norm=norm_l12nn_primal, dual_norm=norm_l12nn_dual, iter_lim=20, verbosity=0) + assert Xnn.dtype == par['dtype'] assert np.any(Xnn < 0) == False @@ -276,6 +343,8 @@ def Aprodfun(A, x, mode): return np.dot(np.conj(A.T), x) return y + np.random.seed(0) + # Create random m-by-n encoding matrix m = par['m'] n = par['n'] diff --git a/spgl1/spgl1.py b/spgl1/spgl1.py index 2294b6c..ecec274 100755 --- a/spgl1/spgl1.py +++ b/spgl1/spgl1.py @@ -45,7 +45,7 @@ def __init__(self, A, nnz_idx, ebar, n): self.shape = (A.shape[0], self.nbar) self.dtype = A.dtype def _matvec(self, x): - y = np.zeros(self.n) + y = np.zeros(self.n, dtype=x.dtype) y[self.nnz_idd] = \ x - (1. / self.nbar) * np.dot(np.dot(np.conj(self.ebar), x), self.ebar) @@ -106,7 +106,7 @@ def _oneprojector_i(b, tau): b = b[idx] csb = np.cumsum(b) - tau - alpha = np.zeros(n+1) + alpha = np.zeros(n + 1) alpha[1:] = csb / (np.arange(n) + 1.0) alphaindex = np.where(alpha[1:] >= b)[0] if alphaindex.any(): @@ -160,11 +160,11 @@ def oneprojector(b, d, tau): Projects b onto the (weighted) one-norm ball of radius tau. If d=1 solves the problem:: - minimize_x ||b-x||_2 st ||x||_1 <= tau. + minimize_x ||b-x||_2 subject to ||x||_1 <= tau. else:: - minimize_x ||b-x||_2 st || Dx ||_1 <= tau. + minimize_x ||b-x||_2 subject to ||Dx||_1 <= tau. Parameters ---------- @@ -881,7 +881,7 @@ def spgl1(A, b, tau=0, sigma=0, x0=None, fid=None, verbosity=0, # Determine initial x and see if problem is complex realx = np.lib.isreal(A).all() and np.lib.isreal(b).all() if x0 is None: - x = np.zeros(n) + x = np.zeros(n, dtype=b.dtype) else: x = np.asarray(x0) @@ -1174,7 +1174,7 @@ def spgl1(A, b, tau=0, sigma=0, x0=None, fid=None, verbosity=0, # LSQR iterations successful. Take the subspace step. if istop != 4: # Push dx back into full space: dx = Z dx. - dx = np.zeros(n) + dx = np.zeros(n, dtype=x.dtype) dx[nnz_idx] = \ dxbar - (1/nebar)*np.dot(np.dot(np.conj(ebar.T), dxbar), dxbar) @@ -1288,7 +1288,7 @@ def spgl1(A, b, tau=0, sigma=0, x0=None, fid=None, verbosity=0, elif stat == EXIT_LEAST_SQUARES: _printf(fid, 'EXIT -- Found a least-squares solution') elif stat == EXIT_LINE_ERROR: - _printf(fid, 'ERROR EXIT -- Linesearch error (%i)' % lnerr) + _printf(fid, 'ERROR EXIT -- Linesearch error (%d)' % lnerr) elif stat == EXIT_SUBOPTIMAL_BP: _printf(fid, 'EXIT -- Found a suboptimal BP solution') elif stat == EXIT_MATVEC_LIMIT: