pub fn difference_matrix(n: usize, order: usize) -> DMatrix<f64>
Compute difference matrix for P-spline penalty.