1
$\begingroup$

I'm doing piece-wise Gaussian linear regression with one breakpoint in raw Python and Numpy and was wondering if it is possible to impose a condition on the weights. Here is the data and sample fit:

enter image description here

I would like the lines to be perpendicular to one another. Mathematically it is simple, the weight of the second line segment should be the negative inverse of the first one: $w_2 = \frac{-1}{w_1}$. However, I am not sure how to incorporate this into the analytical solution for the regression function posterior. Is there a way to do it?

I am currently using the following feature set with the second feature activating at the breakpoint: $$ \phi_x = \begin{bmatrix} 1 \\ x \\ x \end{bmatrix} $$

and following the lectures by Dr. Philipp Hennig (https://www.youtube.com/watch?v=EF1BfKnINw0) with the following analytical formulas:

enter image description here

Here is my code:

    # a is the x breakpoint coordinate, st - x of first data point, corner - x of the breakpoint
    def phi(a, st, corner):
    # ReLU: 
    F = 3
    line_features =  1 * (a - np.array([st[0], corner[0]]).T) * (a > np.array([st[0], corner[0]]).T)
    result = np.ones((line_features.shape[0], line_features.shape[1] + 1))
    result[:, 1:] = line_features
    return result
    
    # first, define the prior
    # number of features
    F = 3
    
    # set parameters of prior on the weights
    mu = np.zeros((F, 1))
    Sigma = 10 * np.eye(F) / F  # p(w)=N(mu,Sigma)
    
    # construct implied prior on f_x
    n = 100  # number of grid-points, for plotting
    
    x = np.linspace(start[0], end[0], n)[:, np.newaxis]  # reshape is needed for phi to work
    
    m = phi(x, start, corner_orig) @ mu
    kxx = phi(x, start, corner_orig) @ Sigma @ phi(x, start, corner_orig).T  # p(f_x)=N(m,k_xx)
    s = multivariate_normal(m.flatten(), kxx + 1e-6 * np.eye(n), size=5).T
    stdpi = np.sqrt(np.diag(kxx))[:, np.newaxis]  # marginal stddev, for plotting
    
    # then, load data from disk
    data = scipy.io.loadmat("nlindata.mat")
    import scipy.io; data = scipy.io.loadmat('nlindata.mat') # use this line to get the nonlinear data.
    X = data["X"]  # inputs
    Y = data["Y"]  # outputs
    sigma = float(data["sigma"])  # measurement noise std-deviation
    
    N = len(X)  # number of data
    
    # evidence: p(Y) = N(Y;M,kXX + sigma**2 * no.eye(N))
    M = phi(X, start, corner_orig) @ mu
    kXX = phi(X, start, corner_orig) @ Sigma @ phi(X, start, corner_orig).T  # p(f_X) = N(M,k_XX)
    
    G = kXX + sigma ** 2 * np.eye(N)
   
    # now, do inference (i.e. construct the posterior)
    # the following in-place decomposition is the most expensive step at O(N^3):
    G = cho_factor(G)
    kxX = phi(x, start, corner_orig) @ Sigma @ phi(X, start, corner_orig).T  # Cov(f_x,f_X) = k_xX
    A = cho_solve(G, kxX.T).T  # pre-compute for re-use (but is only O(N^2))
    
    # # posterior p(f_x|Y) = N(f_x,mpost,vpost)
    mpost = m + A @ (Y - M)  # mean
    vpost = kxx - A @ kxX.T  # covariance
    
    spost = multivariate_normal(mpost.flatten(), vpost  + 1e-6 * np.eye(n), size=5).T  # samples
    stdpo = np.sqrt(np.diag(vpost))[:, np.newaxis]
Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Do you wish to specify the breakpoint or do you need to estimate that, too? You appear to have an interesting version of the problem discussed at stats.stackexchange.com/a/291598/919 wherein you need to impose an additional constraint that $\beta\beta^\prime=1$ (in the notation of that post). Your problem is amenable to the same form of solution (although the algorithm would have to be modified to respect the constraint). $\endgroup$
    – whuber ♦
    Commented Nov 16, 2021 at 17:36
  • $\begingroup$ Thanks for your answer! Yes, I do need a breakpoint estimate (preferably using hierarhical Bayesian inference or other statistical ways that don't rely on optimizations like L2 error reduction) and need to find a way to impose that perpendicularity constraint as well. However, the solution of [link](stats.stackexchange.com/a/291598/919) provides an answer to neither of the issues. There, the breakpoint location is selected manually and the discussion is about finding the best way to fit the data once we know the breakpoint. $\endgroup$
    – Taras Machula
    Commented Nov 17, 2021 at 13:57
  • $\begingroup$ Although your objections apply to the referenced post, my point is that the method there is easily changed to find the breakpoint and impose the constraint. I am puzzled, though, by your requirement that the solution should "not rely on optimizations." That is a crippling restriction. Could you explain why it might be necessary? $\endgroup$
    – whuber ♦
    Commented Nov 17, 2021 at 14:19
  • $\begingroup$ Again, thanks for the answer :) the reason is that I will have to port this to an embedded device with C++ and wanted to avoid additional computational efforts. However, I'm open to suggestions. You mentioned that the solution can be easily adapted to solve my problem, could you please tell how? $\endgroup$
    – Taras Machula
    Commented Nov 17, 2021 at 15:40

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.