I've been looking for an efficient mean-preserving calculation for quadratic interpolation. It doesn't seem to exist. So I tried a bare-bones basic system of linear equations in python and it's numerically unstable even if the input array is completely flat. The model: every bin of data is of width 1. The value of the bin is identical to the integral of a quadratic polynomial over the bin, such that Int(a_i+b_ix+c_ix^2,{dx,x-1/2,x+1/2})=data[x,i, x=i]=a_i+b_ix+c_i*(x^2+1/12). In addition I set values to match between proximate pairs, so that a_i+b_i*(x+1/2)+c_i*(x+1/2)^2=a_{i+1}+b_{i+1}*(x+1/2)+c_{i+1}*(x+1/2)^2 and so on for the slope. The knots of the splines are tied at +-1/2, where the negative end of the first set is left dangling, and I define the boundary conditions on the last set. So for a dataset of length n, I've got 3n equations and 3n variables. np.linalg.solve() can take it! The results are unstable for all input datasets, and immediately unstable for the input dataset [1,2,3]. The system of equations don't get anywhere close to being solved. I expected integro_quadratic_spline_orig to return I after being given the output of integro_quadratic_spline(I). But alas, it comes nowhere near. def integro_quadratic_spline(I): # Number of intervals n = len(I) # Initialize coefficients a = np.zeros(n) b = np.zeros(n) c = np.zeros(n) # Set up the system of equations A = sp.sparse.csr_array((3*n, 3*n)) B = np.zeros(3*n) # Integral conditions #ALL CONDITIONS. First line solve, 2nd line continuity for i in range(0,n-1): A[3*(i),3*(i):3*(i)+6]=[1,(i),(i)**2+1/12,0,0,0] #condition that the integral over the bin matches the original bin value A[3*(i)+1,3*(i):3*(i)+6]=[0,1,2*(i+1/2),0,-1,-2*(i+1/2)] #condition that the slope at the right edge of the current bin matches the slope at the left edge of the next bin A[3*(i)+2,3*(i):3*(i)+6]=[1,i+1/2,(i+1/2)**2,-1,-(i+1/2),-(i+1/2)**2] #condition that the value at the right edge of the current bin matches the value at the left edge of the next bin B[3*(i)]=I #Value of the integral B[3*(i)+1]=0 #matching condition, a'-b'=0 B[3*(i)+2]=0 #matching condition, a-b=0 #below is same as above, but for trailing end. i=i+1 A[3*(i),3*(i):3*(i)+3]=[1,(i),(i)**2+1/12] A[3*(i)+1,3*(i):3*(i)+3]=[0,1,2*(i+1/2)] A[3*(i)+2,3*(i):3*(i)+3]=[1,i+1/2,(i+1/2)**2] B[3*(i)]=I B[3*(i)+1]=0 B[3*(i)+2]=I # Solve the system coeffs = sp.sparse.linalg.spsolve(A, B) print(np.allclose(np.dot(A.toarray(), coeffs), B)) # Extract coefficients for i in range(n): a = coeffs[3*i] b = coeffs[3*i+1] c = coeffs[3*i+2] return a, b, c def integro_quadratic_spline_return(coeffs, mag): original_length = len(coeffs[0]) new_length=original_length*mag newplot=np.empty(new_length) for i in range(original_length): for j in range(int(mag)): newplot[mag*i+j]=coeffs[0]+coeffs[1]*(i+j/mag-1/2)+coeffs[2]*(i+j/mag-1/2)**2 return newplot def integro_quadratic_spline_orig(coeffs): original_length = len(coeffs[0]) newplot=np.empty(original_length) for i in range(original_length): newplot=coeffs[0]+coeffs[1]*(i)+coeffs[2]*((i)**2-12) return newplot Continue reading...
