@@ -231,16 +231,17 @@ def spectraldiff(x, dt, params=None, options=None, high_freq_cutoff=None, even_e
231231
232232
233233def rbfdiff (x , _t , sigma = 1 , lmbd = 0.01 ):
234- """Find smoothed function and derivative estimates by fitting noisy data against a radial-basis-function-filled
235- kernel (naively NxN but made sparse by truncating tiny values). This function uses a Gaussian kernel, thereby
236- assuming the signal is generated from a Gaussian process, effectively the same assumption the Kalman filter makes.
237- It can handle variable step size just as easily as uniform step size.
234+ """Find smoothed function and derivative estimates by fitting noisy data with radial-basis-functions. Naively,
235+ fill a matrix with basis function samples, similar to the implicit inverse problem of spectral methods, but
236+ truncate tiny values to make columns sparse. Each basis function "hill" is topped with a "tower" of height
237+ :code:`lmbd` to reach noisy data samples, and the final smoothed reconstruction is found by razing these and only
238+ keeping the hills.
238239
239240 :param np.array[float] x: data to differentiate
240241 :param float or array[float] _t: This function supports variable step size. This parameter is either the constant
241242 :math:`\\ Delta t` if given as a single float, or data locations if given as an array of same length as :code:`x`.
242- :param float sigma: controls width of radial basis function
243- :param float lmbd: controls strength of bias toward data
243+ :param float sigma: controls width of radial basis functions
244+ :param float lmbd: controls smoothness
244245
245246 :return: tuple[np.array, np.array] of\n
246247 - **x_hat** -- estimated (smoothed) x
@@ -255,10 +256,10 @@ def rbfdiff(x, _t, sigma=1, lmbd=0.01):
255256 # The below does the approximate equivalent of this code, but sparsely in O(N sigma^2), since the rbf falls off rapidly
256257 # t_i, t_j = np.meshgrid(t,t)
257258 # r = t_j - t_i # radius
258- # rbf = np.exp(-(r**2) / (2 * sigma**2)) # radial basis function kernel
259+ # rbf = np.exp(-(r**2) / (2 * sigma**2)) # radial basis function kernel, O(N^2) entries
259260 # drbfdt = -(r / sigma**2) * rbf # derivative of kernel
260261 # rbf_regularized = rbf + lmbd*np.eye(len(t))
261- # alpha = np.linalg.solve(rbf_regularized, x)
262+ # alpha = np.linalg.solve(rbf_regularized, x) # O(N^3)
262263
263264 cutoff = np .sqrt (- 2 * sigma ** 2 * np .log (1e-4 ))
264265 rows , cols , vals , dvals = [], [], [], []
@@ -272,9 +273,9 @@ def rbfdiff(x, _t, sigma=1, lmbd=0.01):
272273 dv = - radius / sigma ** 2 * v
273274 rows .append (n ); cols .append (j ); vals .append (v ); dvals .append (dv )
274275
275- rbf = sparse .csr_matrix ((vals , (rows , cols )), shape = (len (t ), len (t ))) # Build sparse kernels
276+ rbf = sparse .csr_matrix ((vals , (rows , cols )), shape = (len (t ), len (t ))) # Build sparse kernels, O(N sigma) entries
276277 drbfdt = sparse .csr_matrix ((dvals , (rows , cols )), shape = (len (t ), len (t )))
277- rbf_regularized = rbf + lmbd * sparse .eye (len (t ), format = "csr" )
278- alpha = sparse .linalg .spsolve (rbf_regularized , x ) # solve sparse system
278+ rbf_regularized = rbf + lmbd * sparse .eye (len (t ), format = "csr" ) # identity matrix gives a little extra height at the centers
279+ alpha = sparse .linalg .spsolve (rbf_regularized , x ) # solve sparse system targeting the noisy data, O(N sigma^2)
279280
280- return rbf @ alpha , drbfdt @ alpha
281+ return rbf @ alpha , drbfdt @ alpha # find samples of reconstructions using the smooth bases
0 commit comments