function P = penalty_matrix(nf, dcoeff) Inputs: nf: number of features dcoeff: coefficient vector for 0th, 2nd, 3rd derivative and so on Output: P: a [nf]x[nf] symmetric matrix. It contains a weighted sum of different matrices D_i'*D_i. the coefficients of the sum are given by dcoeff. D_0 is eye(nf, nf). D_i is the i-th order differential operator such that D_i*x = diff(x, i).