From bbcf5a339712e5c7ce5cf4e9de6db4ab8bc12f82 Mon Sep 17 00:00:00 2001 From: Philip Loche Date: Fri, 24 Feb 2023 13:35:00 +0100 Subject: [PATCH 1/2] Test an exact solve solver --- src/equisolve/numpy/models/linear_model.py | 42 +++++++++++++++++----- tests/numpy/models/test_linear_model.py | 4 +-- 2 files changed, 34 insertions(+), 12 deletions(-) diff --git a/src/equisolve/numpy/models/linear_model.py b/src/equisolve/numpy/models/linear_model.py index 60ddf67..4e3edf9 100644 --- a/src/equisolve/numpy/models/linear_model.py +++ b/src/equisolve/numpy/models/linear_model.py @@ -16,6 +16,8 @@ from ...utils.metrics import rmse from ..utils import block_to_array, dict_to_tensor_map, tensor_map_to_dict +import scipy.linalg + class Ridge: r"""Linear least squares with l2 regularization for :class:`equistore.Tensormap`'s. @@ -119,18 +121,39 @@ def _validate_params( ) def _numpy_lstsq_solver(self, X, y, sample_weights, alphas, rcond): - # Convert problem with regularization term into an equivalent - # problem without the regularization term + """Solving regularized linear least squares with :class:`numpy.linal.lstsq`. + + The functions converts the problem with regularization term into an equivalent + problem without the regularization term which can be applied to a solver. + """ + num_properties = X.shape[1] - regularization_all = np.hstack((sample_weights, alphas)) + regularization_all = np.hstack((sample_weights[:, 0], alphas[0, :])) regularization_eff = np.diag(np.sqrt(regularization_all)) X_eff = regularization_eff @ np.vstack((X, np.eye(num_properties))) - y_eff = regularization_eff @ np.hstack((y, np.zeros(num_properties))) + y_eff = regularization_eff @ np.hstack((y[:, 0], np.zeros(num_properties))) return np.linalg.lstsq(X_eff, y_eff, rcond=rcond)[0] + def _scipy_solve_solver(self, X, y, sample_weights, alphas): + """Solving regularized linear least squares with :class:`scipy.linal.lstsq`.""" + + X_eff = X / sample_weights + y_eff = y / sample_weights + + X_eff = X.T @ X + np.diag(alphas[0, :]) + y_eff = X.T @ y[:, 0] + + # print(y_eff.shape) + # print(X_eff.shape) + # print(np.diag(alphas[0, :]).shape) + # print(y[:, 0].shape) + + #return scipy.linalg.solve(X_eff, y_eff, assume_a="pos", overwrite_a=True).ravel() + return np.linalg.solve(X_eff, y_eff) + def fit( self, X: TensorMap, @@ -183,23 +206,24 @@ def fit( X_arr = block_to_array(X_block, self.parameter_keys) # y_arr has shape lentgth of n_targets - y_arr = block_to_array(y_block, self.parameter_keys)[:, 0] + y_arr = block_to_array(y_block, self.parameter_keys) # alpha_arr has length of n_properties - alpha_arr = alpha_block.values[0] + alpha_arr = alpha_block.values # Sample weights if sample_weight is not None: sw_block = sample_weight.block(key) # sw_arr has length of n_targets - sw_arr = block_to_array(sw_block, self.parameter_keys)[:, 0] + sw_arr = block_to_array(sw_block, self.parameter_keys) assert ( sw_arr.shape == y_arr.shape ), f"shapes = {sw_arr.shape} and {y_arr.shape}" else: - sw_arr = np.ones((len(y_arr),)) + sw_arr = np.ones((len(y_arr), 1)) - w = self._numpy_lstsq_solver(X_arr, y_arr, sw_arr, alpha_arr, rcond) + #w = self._numpy_lstsq_solver(X_arr, y_arr, sw_arr, alpha_arr, rcond) + w = self._scipy_solve_solver(X_arr, y_arr, sw_arr, alpha_arr) weights_block = TensorBlock( values=w.reshape(1, -1), diff --git a/tests/numpy/models/test_linear_model.py b/tests/numpy/models/test_linear_model.py index 615e13e..1667ea5 100644 --- a/tests/numpy/models/test_linear_model.py +++ b/tests/numpy/models/test_linear_model.py @@ -59,9 +59,7 @@ def to_equistore(self, X_arr=None, y_arr=None, alpha_arr=None, sw_arr=None): def equisolve_solver_from_numpy_arrays(self, X_arr, y_arr, alpha_arr, sw_arr=None): X, y, alpha, sw = self.to_equistore(X_arr, y_arr, alpha_arr, sw_arr) - clf = Ridge( - parameter_keys="values", - ) + clf = Ridge(parameter_keys="values") clf.fit(X=X, y=y, alpha=alpha, sample_weight=sw) return clf From a67fb0edd95ef7c597c9b36bd42f3094bd1fd8e9 Mon Sep 17 00:00:00 2001 From: Philip Loche Date: Wed, 1 Mar 2023 07:57:19 +0000 Subject: [PATCH 2/2] Use scipy solve w --- src/equisolve/numpy/models/linear_model.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/equisolve/numpy/models/linear_model.py b/src/equisolve/numpy/models/linear_model.py index 4e3edf9..d07d4dd 100644 --- a/src/equisolve/numpy/models/linear_model.py +++ b/src/equisolve/numpy/models/linear_model.py @@ -120,7 +120,7 @@ def _validate_params( "properties." ) - def _numpy_lstsq_solver(self, X, y, sample_weights, alphas, rcond): + def _lstsq_solver(self, X, y, sample_weights, alphas, rcond): """Solving regularized linear least squares with :class:`numpy.linal.lstsq`. The functions converts the problem with regularization term into an equivalent @@ -137,7 +137,7 @@ def _numpy_lstsq_solver(self, X, y, sample_weights, alphas, rcond): return np.linalg.lstsq(X_eff, y_eff, rcond=rcond)[0] - def _scipy_solve_solver(self, X, y, sample_weights, alphas): + def _solve_solver(self, X, y, sample_weights, alphas): """Solving regularized linear least squares with :class:`scipy.linal.lstsq`.""" X_eff = X / sample_weights @@ -151,8 +151,7 @@ def _scipy_solve_solver(self, X, y, sample_weights, alphas): # print(np.diag(alphas[0, :]).shape) # print(y[:, 0].shape) - #return scipy.linalg.solve(X_eff, y_eff, assume_a="pos", overwrite_a=True).ravel() - return np.linalg.solve(X_eff, y_eff) + return scipy.linalg.solve(X_eff, y_eff, assume_a="pos", overwrite_a=True).ravel() def fit( self, @@ -161,6 +160,7 @@ def fit( alpha: Union[float, TensorMap] = 1.0, sample_weight: Optional[TensorMap] = None, rcond: float = 1e-13, + solver="solve", ) -> None: """Fit Ridge regression model to each block in X. @@ -221,9 +221,13 @@ def fit( ), f"shapes = {sw_arr.shape} and {y_arr.shape}" else: sw_arr = np.ones((len(y_arr), 1)) - - #w = self._numpy_lstsq_solver(X_arr, y_arr, sw_arr, alpha_arr, rcond) - w = self._scipy_solve_solver(X_arr, y_arr, sw_arr, alpha_arr) + + if solver == "lstsq": + w = self._lstsq_solver(X_arr, y_arr, sw_arr, alpha_arr, rcond) + elif solver == "solve": + w = self._solve_solver(X_arr, y_arr, sw_arr, alpha_arr) + else: + raise ValueError("Unknown solver") weights_block = TensorBlock( values=w.reshape(1, -1),