Source code for equadratures.subspaces

from equadratures.parameter import Parameter
from equadratures.poly import Poly
from equadratures.basis import Basis
from equadratures.scalers import scaler_minmax, scaler_meanvar, scaler_custom
import equadratures.plot as plot
import numpy as np
import scipy
from scipy.linalg import orth, sqrtm
from scipy.spatial import ConvexHull
from scipy.special import comb
from scipy.optimize import linprog
import warnings

[docs]class Subspaces(object): """ This class defines a subspaces object. It can be used for polynomial-based subspace dimension reduction. Parameters ---------- method : str The method to be used for subspace-based dimension reduction. Two options: - ``active-subspace``, which uses ideas in [1] and [2] to compute a dimension-reducing subspace with a global polynomial approximant. Gradients evaluations of the polynomial approximation are used to compute the averaged outer product of the gradient covariance matrix. The polynomial approximation in the original full-space can be provided via ``full_space_poly``. Otherwise, it is fit internally to the data provided via ``sample_points`` and ``sample_outputs``. - ``variable-projection`` [3], where a Gauss-Newton optimisation problem is solved to compute both the polynomial coefficients and its subspace, with the data provided via ``sample_points`` and ``sample_outputs``. full_space_poly : Poly, optional An instance of Poly fitted to the full-space data, to use for the AS computation. sample_points : numpy.ndarray, optional Array with shape (number_of_observations, dimensions) that corresponds to a set of sample points over the parameter space. sample_outputs : numpy.ndarray, optional Array with shape (number_of_observations, 1) that corresponds to model evaluations at the sample points. subspace_dimension : int, optional The dimension of the *active* subspace. param_args : dict, optional Arguments passed to parameters of the AS polynomial. (see :class:`~equadratures.parameter.Parameter`) poly_args : dict , optional Arguments passed to constructing polynomial used for AS computation. (see :class:`~equadratures.poly.Poly`) dr_args : dict, optional Arguments passed to customise the VP optimiser. See documentation for :meth:`~equadratures.subspaces.Subspaces._get_variable_projection` in source. Examples -------- Obtaining a 2D subspace via active subspaces on user data >>> mysubspace = Subspaces(method='active-subspace', sample_points=X, sample_outputs=Y) >>> eigs = mysubspace.get_eigenvalues() >>> W = mysubspace.get_subspace()[:, :2] >>> e = mysubspace.get_eigenvalues() Obtaining a 2D subspace via active subspaces with a Poly object (remember to call set_model() on Poly first) >>> mysubspace = Subspaces(method='active-subspace', full_space_poly=my_poly) >>> eigs = mysubspace.get_eigenvalues() >>> W = mysubspace.get_subspace()[:, :2] >>> e = mysubspace.get_eigenvalues() Obtaining a 2D subspace via variable projection on user data >>> mysubspace = Subspaces(method='variable-projection', sample_points=X, sample_outputs=Y) >>> W = mysubspace.get_subspace()[:, :2] References ---------- 1. Constantine, P., (2015) Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies. SIAM Spotlights. 2. Seshadri, P., Shahpar, S., Constantine, P., Parks, G., Adams, M. (2018) Turbomachinery Active Subspace Performance Maps. Journal of Turbomachinery, 140(4), 041003. `Paper <>`__. 3. Hokanson, J., Constantine, P., (2018) Data-driven Polynomial Ridge Approximation Using Variable Projection. SIAM Journal of Scientific Computing, 40(3), A1566-A1589. `Paper <>`__. """ def __init__(self, method, full_space_poly=None, sample_points=None, sample_outputs=None, subspace_dimension=2, polynomial_degree=2, param_args=None, poly_args=None, dr_args=None): self.full_space_poly = full_space_poly self.sample_points = sample_points self.Y = None # for the zonotope vertices self.sample_outputs = sample_outputs self.method = method self.subspace_dimension = subspace_dimension self.polynomial_degree = polynomial_degree my_poly_args = {'method': 'least-squares', 'solver_args': {}} if poly_args is not None: my_poly_args.update(poly_args) self.poly_args = my_poly_args my_param_args = {'distribution': 'uniform', 'order': self.polynomial_degree, 'lower': -1, 'upper': 1} if param_args is not None: my_param_args.update(param_args) # I suppose we can detect if lower and upper is present to decide between these categories? bounded_distrs = ['analytical', 'beta', 'chebyshev', 'arcsine', 'truncated-gaussian', 'uniform'] unbounded_distrs = ['gaussian', 'normal', 'gumbel', 'logistic', 'students-t', 'studentst'] semi_bounded_distrs = ['chi', 'chi-squared', 'exponential', 'gamma', 'lognormal', 'log-normal', 'pareto', 'rayleigh', 'weibull'] if dr_args is not None: if 'standardize' in dr_args: dr_args['standardise'] = dr_args['standardize'] if self.method.lower() == 'active-subspace' or self.method.lower() == 'active-subspaces': self.method = 'active-subspace' if dr_args is not None: self.standardise = getattr(dr_args, 'standardise', True) else: self.standardise = True if self.full_space_poly is None: # user provided input/output data N, d = self.sample_points.shape if self.standardise: self.data_scaler = scaler_minmax() self.std_sample_points = self.data_scaler.transform(self.sample_points) else: self.std_sample_points = self.sample_points.copy() param = Parameter(**my_param_args) if param_args is not None: if (hasattr(dr_args, 'lower') or hasattr(dr_args, 'upper')) and self.standardise: warnings.warn('Points standardised but parameter range provided. Overriding default ([-1,1])...', UserWarning) myparameters = [param for _ in range(d)] mybasis = Basis("total-order") mypoly = Poly(myparameters, mybasis, sampling_args={'sample-points': self.std_sample_points, 'sample-outputs': self.sample_outputs}, **my_poly_args) mypoly.set_model() self.full_space_poly = mypoly else: # User provided polynomial # Standardise according to distribution specified. Only care about the scaling (not shift) # TODO: user provided callable with parameters? user_params = self.full_space_poly.parameters d = len(user_params) self.sample_points = self.full_space_poly.get_points() if self.standardise: scale_factors = np.zeros(d) centers = np.zeros(d) for dd, p in enumerate(user_params): if in bounded_distrs: scale_factors[dd] = (p.upper - p.lower) / 2.0 centers[dd] = (p.upper + p.lower) / 2.0 elif in unbounded_distrs: scale_factors[dd] = np.sqrt(p.variance) centers[dd] = p.mean else: scale_factors[dd] = np.sqrt(p.variance) centers[dd] = 0.0 self.param_scaler = scaler_custom(centers, scale_factors) self.std_sample_points = self.param_scaler.transform(self.sample_points) else: self.std_sample_points = self.sample_points.copy() if not hasattr(self.full_space_poly, 'coefficients'): raise ValueError('Please call set_model() first on poly.') self.sample_outputs = self.full_space_poly.get_model_evaluations() # TODO: use dr_args for resampling of gradient points as_args = {'grad_points': None} if dr_args is not None: as_args.update(dr_args) self._get_active_subspace(**as_args) elif self.method == 'variable-projection': self.data_scaler = scaler_minmax() self.std_sample_points = self.data_scaler.transform(self.sample_points) if dr_args is not None: vp_args = {'gamma':0.1, 'beta':1e-4, 'tol':1e-7, 'maxiter':1000, 'U0':None, 'verbose':False} vp_args.update(dr_args) self._get_variable_projection(**vp_args) else: self._get_variable_projection()
[docs] def get_subspace_polynomial(self): """ Returns a polynomial defined over the dimension reducing subspace. Returns ------- Poly A Poly object that defines a polynomial over the subspace. The distribution of parameters is assumed to be uniform and the maximum and minimum bounds for each parameter are defined by the maximum and minimum values of the project samples. """ # TODO: Try correlated poly here active_subspace = self._subspace[:, 0:self.subspace_dimension] projected_points =, active_subspace) myparameters = [] for i in range(0, self.subspace_dimension): param = Parameter(distribution='uniform', lower=np.min(projected_points[:, i]), upper=np.max(projected_points[:, i]), order=self.polynomial_degree) myparameters.append(param) mybasis = Basis("total-order") subspacepoly = Poly(myparameters, mybasis, method='least-squares', sampling_args={'sample-points': projected_points, 'sample-outputs': self.sample_outputs}) subspacepoly.set_model() return subspacepoly
[docs] def get_eigenvalues(self): """ Returns the eigenvalues of the dimension reducing subspace. Note: this option is currently only valid for method ``active-subspace``. Returns ------- numpy.ndarray Array of shape (dimensions,) corresponding to the eigenvalues of the above mentioned covariance matrix. """ if self.method == 'active-subspace': return self._eigenvalues else: print('Only the active-subspace method yields eigenvalues.')
[docs] def get_subspace(self): """ Returns the dimension reducing subspace. Returns ------- numpy.ndarray Array of shape (dimensions, dimensions) where the first ``subspace_dimension`` columns contain the dimension reducing subspace, while the remaining columns contain its orthogonal complement. """ return self._subspace
def _get_active_subspace(self, grad_points=None, **kwargs): """ Private method to compute active subspaces. """ if grad_points is None: X = self.full_space_poly.get_points() else: if hasattr(self, 'data_scaler'): X = self.data_scaler.transform(grad_points) else: # Either no standardisation, or user provided poly + param scaling X = grad_points.copy() M, d = X.shape if d != self.sample_points.shape[1]: raise ValueError('In _get_active_subspace: dimensions of gradient evaluation points mismatched with input dimension!') alpha = 2.0 num_grad_lb = alpha * self.subspace_dimension * np.log(d) if M < num_grad_lb: warnings.warn('Number of gradient evaluation points is likely to be insufficient. Consider resampling!', UserWarning) polygrad = self.full_space_poly.get_polyfit_grad(X) if hasattr(self, 'param_scaler'): # Evaluate gradient in transformed coordinate space polygrad = self.param_scaler.div[:, np.newaxis] * polygrad weights = np.ones((M, 1)) / M R = polygrad.transpose() * weights C =, R ) # Compute eigendecomposition! e, W = np.linalg.eigh(C) idx = e.argsort()[::-1] eigs = e[idx] eigVecs = W[:, idx] if hasattr(self, 'data_scaler'): scale_factors = 2.0 / (self.data_scaler.Xmax - self.data_scaler.Xmin) eigVecs = scale_factors[:, np.newaxis] * eigVecs eigVecs = np.linalg.qr(eigVecs)[0] self._subspace = eigVecs self._eigenvalues = eigs def _get_variable_projection(self, gamma=0.1, beta=1e-4, tol=1e-7, maxiter=1000, U0=None, verbose=False): """ Private method to obtain an active subspace in inputs design space via variable projection. Note: It may help to standardize outputs to zero mean and unit variance Parameters ---------- gamma : float, optional Step length reduction factor (0,1). beta : float, optional Armijo tolerance for backtracking line search (0,1). tol : float, optional Tolerance for convergence, measured in the norm of residual over norm of f. maxiter : int, optional Maximum number of optimisation iterations. U0 : numpy.ndarray, optional Initial guess for active subspace. verbose : bool, optional Set to ``True`` for debug messages. """ # NOTE: How do we know these are the best values of gamma and beta? M, m = self.std_sample_points.shape if U0 is None: Z = np.random.randn(m, self.subspace_dimension) U, _ = np.linalg.qr(Z) else: U = orth(U0) y =,U) minmax = np.zeros((2, self.subspace_dimension)) minmax[0, :] = np.amin(y, axis=0) minmax[1, :] = np.amax(y, axis=0) # Construct the affine transformation eta = 2 * np.divide((y - minmax[0,:]), (minmax[1,:]-minmax[0,:])) - 1 # Construct the Vandermonde matrix step 6 V, poly_obj = vandermonde(eta, self.polynomial_degree) V_plus = np.linalg.pinv(V) coeff =, self.sample_outputs) res = self.sample_outputs -,coeff) # R = np.linalg.norm(res) # TODO: convergence criterion?? for iteration in range(0,maxiter): # Construct the Jacobian step 9 J = jacobian_vp(V, V_plus, U, self.sample_outputs, poly_obj, eta, minmax, self.std_sample_points) # Calculate the gradient of Jacobian (step 10) G = np.zeros((m, self.subspace_dimension)) # NOTE: Can be vectorised for i in range(0, M): G += res[i]*J[i, :, :] # conduct the SVD for J_vec vec_J = np.reshape(J, (M, m*self.subspace_dimension)) Y, S, Z = np.linalg.svd(vec_J,full_matrices=False) # step 11 # obtain delta delta =[:,:-self.subspace_dimension**2].T, res) delta =[:-self.subspace_dimension**2]), delta) delta =[:-self.subspace_dimension**2,:].T, delta).reshape(U.shape) # carry out Gauss-Newton step vec_delta=delta.flatten() # step 12 # vectorize G step 13 vec_G = G.flatten() alpha =, vec_delta) norm_G =, vec_G) # check alpha step 14 if alpha >= 0: delta = -G alpha = -norm_G # SVD on delta step 17 Y, S, Z = np.linalg.svd(delta, full_matrices=False) UZ =,Z.T) t = 1 for iter2 in range(0,20): U_new =, np.diag(np.cos(S*t))) +, np.diag(np.sin(S*t)))#step 19 U_new = orth(U_new) # Update the values with the new U matrix y =, U_new) minmax[0,:] = np.amin(y, axis=0) minmax[1,:] = np.amax(y, axis=0) eta = 2 * np.divide((y - minmax[0,:]), (minmax[1,:]-minmax[0,:])) - 1 V_new, poly_obj = vandermonde(eta, self.polynomial_degree) V_plus_new = np.linalg.pinv(V_new) coeff_new =, self.sample_outputs) res_new = self.sample_outputs -,coeff_new) R_new = np.linalg.norm(res_new) if np.linalg.norm(res_new) <= np.linalg.norm(res)+alpha*beta*t or t < 1e-10: # step 21 break t = t * gamma dist_change = subspace_dist(U, U_new) U = U_new V = V_new # coeff = coeff_new V_plus = V_plus_new res = res_new # R = R_new if dist_change < tol: if verbose: print("VP finished with %d iterations" % iteration) break if iteration == maxiter - 1 and verbose: print("VP finished with %d iterations" % iteration) active_subspace = U inactive_subspace = _null_space(active_subspace.T) self._subspace = np.hstack([active_subspace, inactive_subspace])
[docs] def get_zonotope_vertices(self, num_samples=10000, max_count=100000): """ Returns the vertices of the zonotope -- the projection of the high-dimensional space over the computed subspace. Parameters ---------- num_samples : int, optional Number of samples per iteration to check. max_count : int, optional Maximum number of iteration. Returns ------- numpy.ndarray Array of shape (number of vertices, ``subspace_dimension``). Note ---- This routine has been adapted from Paul Constantine's zonotope_vertices() function; see reference below. Constantine, P., Howard, R., Glaws, A., Grey, Z., Diaz, P., Fletcher, L., (2016) Python Active-Subspaces Utility Library. Journal of Open Source Software, 1(5), 79. `Paper <>`__. """ m = self._subspace.shape[0] n = self.subspace_dimension W = self._subspace[:, :n] if n == 1: y0 =, np.sign(W))[0] if y0 < -y0: yl, yu = y0, -y0 xl, xu = np.sign(W), -np.sign(W) else: yl, yu = -y0, y0 xl, xu = -np.sign(W), np.sign(W) Y = np.array([yl, yu]).reshape((2,1)) X = np.vstack((xl.reshape((1,m)), xu.reshape((1,m)))) self.Y = Y return Y else: total_vertices = 0 for i in range(n): total_vertices += comb(m-1,i) total_vertices = int(2*total_vertices) Z = np.random.normal(size=(num_samples, n)) X = get_unique_rows(np.sign(, W.transpose()))) X = get_unique_rows(np.vstack((X, -X))) N = X.shape[0] count = 0 while N < total_vertices: Z = np.random.normal(size=(num_samples, n)) X0 = get_unique_rows(np.sign(, W.transpose()))) X0 = get_unique_rows(np.vstack((X0, -X0))) X = get_unique_rows(np.vstack((X, X0))) N = X.shape[0] count += 1 if count > max_count: break num_vertices = X.shape[0] if total_vertices > num_vertices: print('Warning: {} of {} vertices found.'.format(num_vertices, total_vertices)) Y =, W) self.Y = Y.reshape((num_vertices, n)) return self.Y
[docs] def get_linear_inequalities(self): """ Returns the linear inequalities defining the zonotope vertices, i.e., Ax<=b. Returns ------- tuple Tuple (A,b), containing the numpy.ndarray's A and b; where A is the matrix for setting the linear inequalities, and b is the right-hand-side vector for setting the linear inequalities. """ if self.Y is None: self.Y = self.get_zonotope_vertices() n = self.Y.shape[1] if n == 1: A = np.array([[1],[-1]]) b = np.array([[max(self.Y)],[min(self.Y)]]) return A, b else: convexHull = ConvexHull(self.Y) A = convexHull.equations[:,:n] b = -convexHull.equations[:,n] return A, b
[docs] def get_samples_constraining_active_coordinates(self, inactive_samples, active_coordinates): """ A hit and run type sampling strategy for generating samples at a given coordinate in the active subspace by varying its coordinates along the inactive subspace. Parameters ---------- inactive_samples : int The number of inactive samples required. active_coordiantes : numpy.ndarray The active subspace coordinates. Returns ------- numpy.ndarray Array containing the full-space coordinates. Note ---- This routine has been adapted from Paul Constantine's hit_and_run() function; see reference below. Constantine, P., Howard, R., Glaws, A., Grey, Z., Diaz, P., Fletcher, L., (2016) Python Active-Subspaces Utility Library. Journal of Open Source Software, 1(5), 79. `Paper <>`__. """ y = active_coordinates N = inactive_samples W1 = self._subspace[:, :self.subspace_dimension] W2 = self._subspace[:, self.subspace_dimension:] m, n = W1.shape s =, y).reshape((m, 1)) normW2 = np.sqrt(np.sum(np.power(W2, 2), axis=1)).reshape((m, 1)) A = np.hstack((np.vstack((W2, -W2.copy())), np.vstack((normW2, normW2.copy())))) b = np.vstack((1 - s, 1 + s)).reshape((2 * m, 1)) c = np.zeros((m - n + 1, 1)) c[-1] = -1.0 # print() zc = linear_program_ineq(c, -A, -b) z0 = zc[:-1].reshape((m - n, 1)) # define the polytope A >= b s =, y).reshape((m, 1)) A = np.vstack((W2, -W2)) b = np.vstack((-1 - s, -1 + s)).reshape((2 * m, 1)) # tolerance ztol = 1e-6 eps0 = ztol / 4.0 Z = np.zeros((N, m - n)) for i in range(N): # random direction bad_dir = True count, maxcount = 0, 50 while bad_dir: d = np.random.normal(size=(m - n, 1)) bad_dir = np.any(, z0 + eps0 * d) <= b) count += 1 if count >= maxcount: Z[i:, :] = np.tile(z0, (1, N - i)).transpose() yz = np.vstack([np.repeat(y[:, np.newaxis], N, axis=1), Z.T]) return, yz).T # find constraints that impose lower and upper bounds on eps f, g = b -, z0),, d) # find an upper bound on the step min_ind = np.logical_and(g <= 0, f < -np.sqrt(np.finfo(np.float).eps)) eps_max = np.amin(f[min_ind] / g[min_ind]) # find a lower bound on the step max_ind = np.logical_and(g > 0, f < -np.sqrt(np.finfo(np.float).eps)) eps_min = np.amax(f[max_ind] / g[max_ind]) # randomly sample eps eps1 = np.random.uniform(eps_min, eps_max) # take a step along d z1 = z0 + eps1 * d Z[i, :] = z1.reshape((m - n,)) # update temp var z0 = z1.copy() yz = np.vstack([np.repeat(y[:, np.newaxis], N, axis=1), Z.T]) return, yz).T
[docs] def plot_sufficient_summary(self, ax=None, X_test=None, y_test=None, show=True, poly=True, uncertainty=False, legend=False, scatter_kwargs={}, plot_kwargs={}): """ Generates a sufficient summary plot for 1D or 2D polynomial ridge approximations. See :meth:`~equadratures.plot.plot_sufficient_summary` for full description. """ return plot.plot_sufficient_summary(self, ax, X_test, y_test, show, poly, uncertainty, legend, scatter_kwargs, plot_kwargs)
[docs] def plot_2D_contour_zonotope(self, mysubspace, minmax=[- 3.5, 3.5], grid_pts=180, show=True, ax=None): """ Generates a 2D contour plot of the polynomial ridge approximation. See :meth:`~equadratures.plot.plot_2D_contour_zonotope` for full description. """ return plot.plot_2D_contour_zonotope(self,minmax,grid_pts,show,ax)
[docs] def plot_samples_from_second_subspace_over_first(self, mysubspace_2, axs=None, no_of_samples=500, minmax=[- 3.5, 3.5], grid_pts=180, show=True): """ Generates a zonotope plot where samples from the second subspace are projected over the first. See :meth:`~equadratures.plot.plot_samples_from_second_subspace_over_first` for full description. """ return plot.plot_samples_from_second_subspace_over_first(self,mysubspace_2, axs, no_of_samples, minmax, grid_pts, show)
def vandermonde(eta, p): # TODO: Try using a "correlated" basis here? _, n = eta.shape listing = [] for i in range(0, n): listing.append(p) Object=Basis('total-order',listing) # Establish n Parameter objects params = [] P = Parameter(order=p, lower=-1, upper=1, distribution='uniform') for i in range(0, n): params.append(P) # Use the params list to establish the Poly object poly_obj = Poly(params, Object, method='least-squares') V = poly_obj.get_poly(eta) V = V.T return V, poly_obj def vector_AS(list_of_polys, R = None, alpha=None, k=None, samples=None, bootstrap=False, bs_trials = 50 , J = None, save_path = None): # Find AS directions to vector val func # analogous to computeActiveSubspace # Since we are dealing with *one* vector val func we should have just one input space # Take the first of the polys. poly = list_of_polys[0] if samples is None: d = poly.dimensions if alpha is None: alpha = 4 if k is None or k > d: k = d M = int(alpha * k * np.log(d)) X = np.zeros((M, d)) for j in range(0, d): X[:, j] = np.reshape(poly.parameters[j].getSamples(M), M) else: X = samples M, d = X.shape n = len(list_of_polys) # number of outputs if R is None: R = np.eye(n) elif len(R.shape) == 1: R = np.diag(R) if J is None: J = jacobian_vec(list_of_polys,X) if not(save_path is None):,J) J_new = np.matmul(sqrtm(R), np.transpose(J,[2,0,1])) JtJ = np.matmul(np.transpose(J_new,[0,2,1]), J_new) H = np.mean(JtJ,axis=0) # Compute P_r by solving generalized eigenvalue problem... # Assume sigma = identity for now e, W = np.linalg.eigh(H) eigs = np.flipud(e) eigVecs = np.fliplr(W) if bootstrap: all_bs_eigs = np.zeros((bs_trials, d)) all_bs_W = [] for t in range(bs_trials): print("Starting bootstrap trial %d"%t) bs_samples = X[np.random.randint(0,M,size=M), :] J_bs = jacobian_vec(list_of_polys, bs_samples) J_new_bs = np.matmul(sqrtm(R), np.transpose(J_bs,[2,0,1])) JtJ_bs = np.matmul(np.transpose(J_new_bs, [0, 2, 1]), J_new_bs) H_bs = np.mean(JtJ_bs, axis=0) # Compute P_r by solving generalized eigenvalue problem... # Assume sigma = identity for now e_bs, W_bs = np.linalg.eigh(H_bs) all_bs_eigs[t,:] = np.flipud(e_bs) eigVecs_bs = np.fliplr(W_bs) all_bs_W.append(eigVecs_bs) eigs_bs_lower = np.min(all_bs_eigs, axis = 0) eigs_bs_upper = np.max(all_bs_eigs, axis = 0) return eigs,eigVecs,eigs_bs_lower,eigs_bs_upper, all_bs_W else: return eigs,eigVecs def jacobian_vp(V, V_plus, U, f, Polybasis, eta, minmax, X): M, N = V.shape m, n = U.shape Gradient = Polybasis.get_poly_grad(eta) sub = (minmax[1,:]-minmax[0,:]).T vectord = np.reshape(2.0/sub,(n,1)) # Initialize the tensor J = np.zeros((M, m, n)) # Obtain the derivative of this tensor dV = np.zeros((m, n, M, N)) for l in range(0, n): for j in range(0, N): current = Gradient[l].T if n == 1: current = Gradient.T dV[:,l,:,j] = np.asscalar(vectord[l])*(X.T*current[:,j]) # Get the P matrix P = np.identity(M)-np.matmul(V,V_plus) V_minus = scipy.linalg.pinv(V) # Calculate entries for the tensor for j in range(0,m): for k in range(0,n): temp1 = np.linalg.multi_dot([P,dV[j,k,:,:],V_minus]) J[:, j, k]=(-np.matmul((temp1+temp1.T),f)).reshape((M,)) # Eqn 15 return J def jacobian_vec(list_of_poly, X): m = len(list_of_poly) [N, d] = X.shape J = np.zeros((m, d, N)) for p in range(len(list_of_poly)): J[p, :, :] = list_of_poly[p].get_polyfit_grad(X) return J def subspace_dist(U, V): if len(U.shape) == 1: return np.linalg.norm(np.outer(U, U) - np.outer(V, V), ord=2) else: return np.linalg.norm(, U.T) -, V.T), ord=2) def linear_program_ineq(c, A, b): c = c.reshape((c.size,)) b = b.reshape((b.size,)) # make unbounded bounds bounds = [] for i in range(c.size): bounds.append((None, None)) A_ub, b_ub = -A, -b res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, options={"disp": False}, method='simplex') if res.success: return res.x.reshape((c.size, 1)) else: np.savez('bad_scipy_lp_ineq_{:010d}'.format(np.random.randint(int(1e9))), c=c, A=A, b=b, res=res) raise Exception('Scipy did not solve the LP. Blame Scipy.') def get_unique_rows(X0): X1 = X0.view(np.dtype((np.void, X0.dtype.itemsize * X0.shape[1]))) return np.unique(X1).view(X0.dtype).reshape(-1, X0.shape[1]) def _null_space(A, rcond=None): """ null space method adapted from scipy. """ u, s, vh = scipy.linalg.svd(A, full_matrices=True) M, N = u.shape[0], vh.shape[1] if rcond is None: rcond = np.finfo(s.dtype).eps * max(M, N) tol = np.amax(s) * rcond num = np.sum(s > tol, dtype=int) Q = vh[num:,:].T.conj() return Q