diff --git a/Reproducibility/Lipschitz/fig_ridge.py b/Reproducibility/Lipschitz/fig_ridge.py
new file mode 100644
index 0000000..a6f5434
--- /dev/null
+++ b/Reproducibility/Lipschitz/fig_ridge.py
@@ -0,0 +1,46 @@
+import numpy as np
+import psdr, psdr.demos
+
+np.random.seed(0)
+
+#fun = psdr.demos.OTLCircuit()
+fun = psdr.demos.Borehole()
+
+X = fun.domain.sample(1e2)
+grads = fun.grad(X)
+
+lip = psdr.LipschitzMatrix(verbose = True)
+lip.fit(grads = grads)
+
+
+U = lip.U[:,0:1]
+print(U)
+
+# setup an experimental design
+X = psdr.minimax_design_1d(fun.domain, 20, L = U.T)
+X2 = []
+for x in X:
+	dom = fun.domain.add_constraints(A_eq = U.T, b_eq = U.T @ x)
+	X2.append(dom.sample_boundary(50))
+
+X2 = np.vstack(X2)
+
+fX2 = fun(X2)
+
+pra = psdr.PolynomialRidgeApproximation(9, 1, norm = np.inf)
+
+pra.fit_fixed_subspace(X2, fX2, U)
+
+
+if True:
+	X3 = psdr.minimax_design_1d(fun.domain, 500, L = U.T)
+	yapprox = pra(X3)
+	y = fun(X3)
+	
+	import matplotlib.pyplot as plt
+	fig, ax = plt.subplots()
+	ax.plot((U.T @ X2.T).flatten(), fX2, 'k.')
+	ax.plot((U.T @ X3.T).flatten(), yapprox, 'r-')
+	
+	plt.show()
+
diff --git a/psdr/domains/euclidean.py b/psdr/domains/euclidean.py
index ee97707..9372c39 100644
--- a/psdr/domains/euclidean.py
+++ b/psdr/domains/euclidean.py
@@ -662,6 +662,27 @@ def sample_grid(self, n):
 		I = self.isinside(Xgrid)
 		return Xgrid[I]	
 
+
+	def sample_boundary(self, draw):
+		r""" Sample points on the boundary of the domain
+
+		Parameters
+		----------
+		draw : int
+			Number of points to sample
+		"""
+		draw = int(draw)
+	
+		dirs = [self.random_direction(self.center) for i in range(draw)]
+		X = np.array([self.corner(d) for d in dirs])
+	
+		if draw == 1:
+			return X.flatten()
+		
+		return X
+	
+
+
 	def random_direction(self, x):
 		r""" Returns a random direction that can be moved and still remain in the domain
 
diff --git a/psdr/polyridge.py b/psdr/polyridge.py
index e32bec6..ea99544 100644
--- a/psdr/polyridge.py
+++ b/psdr/polyridge.py
@@ -117,7 +117,7 @@ def inf_norm_fit(A, b):
 	with warnings.catch_warnings():
 		warnings.simplefilter('ignore', PendingDeprecationWarning)
 		x = cp.Variable(A.shape[1])
-		obj = cp.norm_inf(x.__rmatmul__(A) - b)
+		obj = cp.norm_inf(A @ x - b.flatten())
 		problem = cp.Problem(cp.Minimize(obj))
 		problem.solve(solver = 'ECOS')
 		return x.value
@@ -295,6 +295,16 @@ def __len__(self):
 	def __str__(self):
 		return "<PolynomialRidgeApproximation degree %d, subspace dimension %d>" % (self.degree, self.subspace_dimension)
 
+
+	def fit_fixed_subspace(self, X, fX, U):
+		r"""
+
+		"""
+		assert U.shape[0] == X.shape[1], "U has %d rows, expected %d based on X" % (U.shape[0], X.shape[1])
+		assert U.shape[1] == self.subspace_dimension, "U has %d columns; expected %d" % (U.shape[1], self.subspace_dimension)
+		self._finish(X, fX, U)
+		
+
 	def fit(self, X, fX, U0 = None):
 		r""" Given samples, fit the polynomial ridge approximation.