Updated numerical experiments

Reproducibility/Stabilized_SK/fig_pbeam.py

@@ -29,7 +29,7 @@
 Xtrain = X[I,:]
 ytrain = y[I]
 
-paaa = ParametricAAARationalApproximation(tol = 1e-7, maxiter = 6)
+paaa = ParametricAAARationalApproximation(tol = 1e-7, maxiter = 5)
 paaa.fit(Xtrain, ytrain)
 
 num_degree = paaa.num_degree

Reproducibility/Stabilized_SK/fig_penzl.py

@@ -3,7 +3,7 @@
 import tqdm
 import scipy.io
 from polyrat import *
-
+from pgf import PGF
 
 def penzl2(X):
 
@@ -21,30 +21,14 @@ def penzl2(X):
 		p = x[1]
 		A1 = np.array([[-1,p],[-p,-1]])
 		A = block_diag(A1, A2, A3, A4)
-		H[i] = c.T @ np.linalg.solve(z*np.eye(1006) - A, b)
+		#H[i] = c.T @ np.linalg.solve(z*np.eye(1006) - A, b)
+		H[i] += c[0:2].T @ np.linalg.solve(z*np.eye(2) - A1, b[0:2])
+		H[i] += c[2:4].T @ np.linalg.solve(z*np.eye(2) - A2, b[2:4])
+		H[i] += c[4:6].T @ np.linalg.solve(z*np.eye(2) - A3, b[4:6])
+		H[i] += c[6:].T @ (1./(z - np.diag(A4))).reshape(-1,1) 
 
 	return H
 
-def penzl3(X):
-
-	b = np.ones(1006)
-	b[0:6] *= 10
-	c = np.copy(b)
-
-	A4 = -np.diag(np.arange(1, 1001))
-
-	H = np.zeros(len(X), dtype = np.complex)
-	for i, x in enumerate(X):
-		z = x[0]
-		p1 = x[1]
-		p2 = x[2]
-		A1 = np.array([[-1,p1],[-p1,-1]])
-		A2 = np.array([[-1, p2],[-p2,-1]])
-		A3 = np.array([[-1,2*p2], [-2*p2,-1]])
-		A = block_diag(A1, A2, A3, A4)
-		H[i] = c.T @ np.linalg.solve(z*np.eye(1006) - A, b)
-
-	return H
 
 try:
 	dat = scipy.io.loadmat('penzl2.dat')	
@@ -72,22 +56,62 @@ def penzl3(X):
 paaa = ParametricAAARationalApproximation(maxiter= 9)
 paaa.fit(X, y)
 
+
 num_degree = paaa.num_degree
 denom_degree = paaa.denom_degree
-ssk = SKRationalApproximation(num_degree, denom_degree, refine = False, maxiter = 5)
+ssk = SKRationalApproximation(num_degree, denom_degree, refine = False, maxiter = 20)
 ssk.fit(X, y)
 
+print("error on training set P-AAA", np.linalg.norm(paaa(X) - y))
+print("error on training set S-SK", np.linalg.norm(ssk(X) - y))
 
 import matplotlib.pyplot as plt
 fig, ax = plt.subplots(1, 2)
 
 err = np.abs(paaa(Xtest) - ytest).reshape(300,1000)
+
+pgf = PGF()
+pgf.add('z', np.log10(Xtest[:,0].imag))
+pgf.add('s', Xtest[:,1].real)
+pgf.add('err', np.log10(err.flatten()))
+pgf.write('data/fig_penzl_paaa.dat')
+
+# small
+pgf = PGF()
+pgf.add('z', np.log10(X[:,0].imag))
+pgf.add('s', X[:,1].real)
+pgf.add('err', np.log10(np.abs(paaa(X) - y)+1e-50))
+pgf.write('data/fig_penzl_paaa_small.dat')
+
+
 ax[0].set_yscale('log')
-cs = ax[0].contourf(np.abs(S), np.abs(Z), np.log10(err), levels = np.arange(-10,1))
+cs = ax[0].contourf(np.abs(S), np.abs(Z), np.log10(err), levels = np.arange(-8,0))
+
+
+Xi = paaa.interpolation_points
+pgf = PGF()
+pgf.add('z', np.log10(Xi[:,0].imag))
+pgf.add('s', Xi[:,1].real)
+pgf.write('data/fig_penzl_paaa_interp.dat')
+
+ax[0].plot(np.abs(Xi[:,1]), np.abs(Xi[:,0]), 'r.')
 
 err = np.abs(ssk(Xtest) - ytest).reshape(300,1000)
+
+pgf = PGF()
+pgf.add('z', np.log10(Xtest[:,0].imag))
+pgf.add('s', Xtest[:,1].real)
+pgf.add('err', np.log10(err.flatten()))
+pgf.write('data/fig_penzl_ssk.dat')
+
+pgf = PGF()
+pgf.add('z', np.log10(X[:,0].imag))
+pgf.add('s', X[:,1].real)
+pgf.add('err', np.log10(np.abs(ssk(X) - y)+1e-50))
+pgf.write('data/fig_penzl_ssk_small.dat')
+
 ax[1].set_yscale('log')
-ax[1].contourf(np.abs(S), np.abs(Z), np.log10(err), levels = np.arange(-10,1))
+ax[1].contourf(np.abs(S), np.abs(Z), np.log10(err), levels = np.arange(-8,0))
 fig.colorbar(cs)
 plt.show()
 

Reproducibility/Stabilized_SK/fig_penzl3.py

@@ -23,7 +23,11 @@ def penzl3(X):
 		A2 = np.array([[-1, p2],[-p2,-1]])
 		A3 = np.array([[-1,2*p2], [-2*p2,-1]])
 		A = block_diag(A1, A2, A3, A4)
-		H[i] = c.T @ np.linalg.solve(z*np.eye(1006) - A, b)
+		#H[i] = c.T @ np.linalg.solve(z*np.eye(1006) - A, b)
+		H[i] += c[0:2].T @ np.linalg.solve(z*np.eye(2) - A1, b[0:2])
+		H[i] += c[2:4].T @ np.linalg.solve(z*np.eye(2) - A2, b[2:4])
+		H[i] += c[4:6].T @ np.linalg.solve(z*np.eye(2) - A3, b[4:6])
+		H[i] += c[6:].T @ (1./(z - np.diag(A4))).reshape(-1,1) 
 
 	return H
 
@@ -42,12 +46,34 @@ def penzl3(X):
 	scipy.io.savemat('penzl3.dat', {'X':X, 'y':y})
 
 
-paaa = ParametricAAARationalApproximation()
-paaa.fit(X, y)
+d1 = []
+d2 = []
+d3 = []
+paaa_err = []
+ssk_err = []
 
-num_degree = paaa.num_degree
-denom_degree = paaa.denom_degree
-ssk = SKRationalApproximation(num_degree, denom_degree, refine = False, maxiter = 5)
-ssk.fit(X, y)
+for maxiter in range(7, 15):
+	paaa = ParametricAAARationalApproximation(maxiter = maxiter)
+	paaa.fit(X, y)
+	paaa_err.append(np.linalg.norm(paaa(X) - y))
 
+	num_degree = paaa.num_degree
+	denom_degree = paaa.denom_degree
 
+	d1.append(num_degree[0])
+	d2.append(num_degree[1])
+	d3.append(num_degree[2])
+
+
+	ssk = SKRationalApproximation(num_degree, denom_degree, refine = False, maxiter = 20)
+	ssk.fit(X, y)
+	ssk_err.append(np.linalg.norm(ssk(X) - y))
+
+	pgf = PGF()
+	pgf.add('d1', d1)
+	pgf.add('d2', d2)
+	pgf.add('d3', d3)
+	pgf.add('paaa_err', paaa_err)
+	pgf.add('ssk_err', ssk_err)
+	pgf.write('data/fig_penzl3.dat')
+	break