Wolfram Cloud Document

Abstract

The contributions here are to:

◼

show numerical equivalence of three methods of linear regression: MAP (Maximum A-Posteriori), KAL (Kalman filtering/folding), and RLS (Recurrent Least Squares) when they are regularized by Bayesian a-priori beliefs, that is, protected from overfitting

◼

remove risky matrix inversion from MAP, KAL, and RLS

◼

remove expensive storage of observation vectors and matrices from MAP by rewriting MAP in recurrent form

◼

show numerical equivalence of KAL to recurrent MAP when MAP’s regularization hyperparameters

α=1

and

β=1

are swapped and inverted

◼

conjecture formal equivalence and support the conjecture algebraically through order three

◼

sketch application of MAP, KAL, and RLS to Finite-Difference Policy Gradient in Reinforcement Learning

◼

apply all methods to a non-linear Padé approximant

Introduction

We investigate four popular methods for linear regression, i.e., parameter estimation, specifically with respect to the listed properties:

◼

MLE — Maximum Likelihood Estimation

◼

Is exposed to overfitting.

◼

Stores and inverts large matrices.

◼

Has numerical hazards of matrix inversion.

◼

MAP — Maximum A-Posteriori

◼

Addresses overfitting via regularization hyperparameters.

◼

The standard algorithm stores and inverts large matrices.

◼

When MAP’s hyperparameters are swapped and inverted, MAP’s recurrent form resembles KAL.

◼

KAL — Kalman filtering

◼

Addresses overfitting via frank priors: estimates and covariances.

◼

Runs in constant memory. The recurrent form is standard.

◼

Dual co-Kalman and co-RLS address Finite-Difference Policy Gradient in Reinforcement Learning.

◼

RLS — Recurrent Least Squares

◼

Matches KAL when renormalized by observation covariance.

◼

Has apparently lower numerical risk than KAL.

◼

Avoids subtraction in the covariance update.

We show that recurrent MAP, KAL, and RLS are numerically equivalent for appropriate choices of a-priori covariances. We conjecture that they’re mathematically equivalent and present analytical evidence for the conjecture.

We note in passing that KAL and RLS admit full covariance (information) matrices, whilst MAP admits only standard deviations. MAP thus assumes that parameters and observations are non-dimensional and commensurate in scale. This means that practical application of MAP may require extra preconditioning steps. We do not dwell on this point in this article.

We show how to apply MAP, KAL, and RLS to an important practical problem: Finite-Difference Policy Gradient in Reinforcement Learning.

We close with a new example: Padé approximants. Linearization opens them to KAL and RLS filtering.

Overfitting and Regularization

Over-fitting is a hazard of linear regression. Linear models, as their number of parameters nears or exceeds the number of data points, follow the data too well — “wiggle too much.” Such wiggling limits the utility of predictions. Between data points and outside data ranges, predictions become unreliable. We show visual demonstrations below after setting things up.

Regularization is the usual fix for over-fitting. Bayesian MAP regularizes via a-priori belief hyperparameters

α=1

and

β=1

— reciprocal variances of the unknown parameters

and the noise on observations

. MAP does not include prior estimates for the parameters ξ themselves.

Kalman filtering (KAL) and recurrent least squares (RLS) are structurally superior because they incorporate priors as extra observations with covariances. KAL and RLS can be stopped and restarted at any point in an infinite stream of observations, and do not need tuning.

RLS, renormalized with observation covariance, reproduces KAL and MAP within numerical errors. Renormalized RLS can have better numerical stability than KAL, trading off risk of ill-conditioning against risk of catastrophic cancelation.

We show below, algebraically through order three and numerically for orders higher than that, that MAP produces the same estimates as do KAL and renormalized RLS when

and

are swapped and inverted. After swapping and inverting α and β, the MAP equations resemble the equations of Kalman filtering, at least for the estimate. In post-processing, we un-swap and re-invert α and β to recover covariances from recurrent MAP.

MAP, KAL, and RLS all regularize, just in different ways.

KAL and RLS scale better than non-recurrent, full-storage MAP, which is a Bayesian modification of non-recurrent, full-storage MLE. There is no obvious way to convert MLE and MAP into recurrences, let alone to extend them to lazy sequences or asynchronous streams. KAL and RLS, however, are overtly recurrent. They are foldable: independent of the manner and means of data delivery. They avoid matrix storage, inversions, and multiplication.

Motivating Example

We begin with MLE on a familiar school example: estimate the slope

and intercept

of the best-fit line function

z(x)=mx+b

to noisy data

(

)

Think of

z(x)=mx+b

, more generally; as a linear combination of potentially non-linear basis functions of

with coefficients

and

. In the expression

z(x)=mx+b

, there are two basis functions. The first is

(x)=x

, the identity function of

. The linear combination

(x)=mx

happens, accidentally, to be linear in the independent variable

. Do not be distracted by that linearity; it is not the linearity we’re looking for. The linearity of

(x)=mx

in the parameter

is the interesting linearity. The second basis function,

(x)=1

, is the constant 1, non-linear in

. The linear combination

(x)=b×1=b

is linear in the parameter

. Thus the entire linear combination

z(x)=m

(x)+b

(x)=mx+b

is linear in

and

, and all methods of estimating

and

are called linear regression.

The idea of linear combinations of non-linear functions extends to sophisticated models. A Fourier series

∞

∑

k=0

is linear in the amplitudes

, but non-linear in the basis functions

. A linear combination of non-linear Bessel’s functions approximates a vibrating drumhead. A linear combination of non-linear spherical harmonics approximates a lumpy Earth or a Hydrogen atom. These examples emphasize that linear regression methods like MLE, MAP, KAL, and RLS pertain to sophisticated models, linear in the parameters but non-linear in the independent variables.

MLE, Four Ways

We show numerical equivalence of four ways to do MLE:

(WBI) with Wolfram’s built-in functions

(CNE) directly through the classic normal equations

(LPI) via the Moore-Penrose left pseudoinverse

(SLS) sidestepping the (hazardous?) inverse by solving a linear system

◼

Problem Statement

Find best-fitting parameters

(slope) and

(intercept), where

z(x)=mx+b

, given known, noisy data

(

,…,

)

and abscissae

(

,…,

)

. Best-fitting means minimizing mean squared error.

Write this system as a matrix equation:

◼

(observations; known, concrete, numerical; scalars in this example, but generalizable to vectors)

◼

(

partials

, Jacobian of the observation model

z(x)=mx+b

with respect to the unknown parameters

Ξ=(m,b)

; rows containing values of basis functions



i∈[1..k]

,1

; known, concrete, numerical)

◼

Ξ=(m,b)

(model, aka state; unknown abstract, symbolic parameters to be estimated).

rows of data

and

rows of partials

come in matched pairs.

N×1

⋮

x 1	1
x 2	1
⋮	⋮
x N	1

unknown

+noise

def

N×2

2×1

+samplesofNormalDistribution[0,

]

(

)

◼

Ground Truth

Fake some data by (1) sampling a line specified by ground truth parameters

and

, then (2) adding Gaussian observation noise. Run the faked data through the four estimation procedures WBI, CNE, LPI, SLS, and see how close the estimated

estimated

and

estimated

come to ground truth. This model is of order

M=2

. There are

N=119

observed scalar data.

is a column vector of dimension 119, a number picked arbitrarily.

In real-world applications, we rarely have ground truth. Its purpose here is to calibrate the four estimation procedures WBI, CNE, LPI, SLS.

Note and remember the ground-truth values

=0.5,

=-1/3

In[]:=

ClearAll[groundTruth,m,b];groundTruth={m,b}={0.5,-1./3.};

◼

Partials

The partials inhabit an

N×M

matrix

, an order-

column vector of order-

row-co-vectors. Each co-vector is the Jacobian gradient 1-form of

A·Ξ

with respect to

, the vector of unknown parameters, evaluated at specific values of

from the data. Leave partials as a global variable for the following demonstrations.

◼

A major defect of MLE is that

and

are large, potentially infinite, in row count. This defect is solved by recurrence, a major theme of this article.

In[]:=

ClearAll[nData,min,max];

nData=119

;min=-1.;max=3.;ClearAll[partials];

partials

=Array[{#,1.0}&,nData,{min,max}];Short[partials,3]

Out[]//Short=

{{-1.,1.},{-0.966102,1.},{-0.932203,1.},{-0.898305,1.},{-0.864407,1.},{-0.830508,1.},108,{2.86441,1.},{2.89831,1.},{2.9322,1.},{2.9661,1.},{3.,1.}}

◼

Faked Observations Ζ

Define a faking function and a global variable, data, for later use. Assume that the observations have unbiased (zero-mean) Gaussian white noise with variance

=0.65

, a convenient but arbitrary value.

In[]:=

ClearAll[fake];fake[n_,σ_,A_,{m_,b_}]:=Table[RandomVariate[NormalDistribution[0,σ]]+A〚i〛.{m,b},{i,n}];

In[]:=

ClearAll[data,noiseσ];noiseσ=0.65;data=fake[nData,noiseσ,partials,groundTruth];Short[data,3]

Out[]//Short=

{-2.39451,-0.461069,0.247936,-0.694994,-0.26789,-1.98306,107,0.349891,0.848294,0.211755,0.918336,1.93958,0.569446}

◼

Wolfram Built-In

The Wolfram built-in LinearModelFit computes an MLE for

Ξ=



. The results in the cell below should approximate ground truth

=0.5,

=-1/3

reasonably well (the results are randomized every time the notebook is evaluated).

In[]:=

ClearAll[model];model=

LinearModelFit

[{partials〚All,1〛,data},x,x];Normal[model]

Out[]=

-0.323941+0.505112x

Here is everything Wolfram has to say about this MLE (it’s a lot; use the notebook’s interactive buttons to expand and contract the display).

In[]:=

Association[(#model[#])&/@model["Properties"]]

Out[]=

AdjustedRSquared0.405362,AIC260.525,

⋯61⋯

,VarianceInflationFactors{0.,1.}

Full expression not available

(

original memory size:

185.9 kB)

The most important attribute of the model is its covariance matrix. We come back to it below.

In[]:=

model["CovarianceMatrix"]//MatrixForm

Out[]//MatrixForm=



0.00738074	-0.00313283
-0.00313283	0.00313283



The following plot shows that Wolfram does a god job of estimating the parameters

and

. We have 119 data and two parameters to estimate, so overfitting is not a hazard here because the number of data greatly exceeds the number of parameters. We’ll show overfitting soon enough.

◼

Normal Equations

That gives numerically the same answer as Wolfram’s built-in, just in the opposite order:

◼

Moore-Penrose Left PseudoInverse

◼

Avoiding Inverse

In general, we advise replacing any occurrence of Inverse[A].B with LinearSolve[A,B]. See https://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/.

Do Not Use the Normal Equations

Recurrence

Recurrent Least Squares (RLS)

Derive, illustrate, and check a recurrent form of MLE. Regularize it (avoid overfitting) via Bayesian priors.

where

Derivation Sketch

where

Numerical Demonstration

Structural Notes

The highlighted mappings of List over the data and partials lists convert them from 1-dimensional puns of column vectors into 2-D column vectors.

Memory and Time Efficiency

Because the observations and partial are the third argument of Fold, they need not be stored. They might as well be infinite and / or asynchronous streams. See https://community.wolfram.com/groups/-/m/t/3078712 and https://community.wolfram.com/groups/-/m/t/3074231.

Check the a-priori

Covariance of the Estimate

Check that this covariance approximately equals that deduced from Π, the final value of Λ:

This is the same covariance matrix that Wolfram’s LinearModel reports. This second check of our recurrence succeeds.

Covariance of the Prediction

Interim Conclusions

We reduce numerical risk by solving a linear system instead of inverting a matrix.

We still have work to do with renormalization and observation covariances.

Regularization and Renormalization

Chris Bishop’s Pattern Recognition and Machine Learning has an extended example of fitting polynomials starting in section 1.1. The higher the order of the polynomial, the more MLE over-fits. Bishop presents MAP regularization as a cure for this over-fitting. KAL and RLS regularize by construction. In this section, we relate KAL and RLS regularization via frank priors to Bishop’s MAP regularization via hyperparameters.

◼

Reproducing Bishop’s Example

Fit polynomials of various order to a sin wave.

Bishop’s Training Set

◼

Wolfram’s built-in NormalDistribution takes the standard deviation as its second argument, not the variance. Bishop’s notation for normal distribution takes variance as second argument, so beware.

Take a sample of the outputs and assign it the name bts for bishopTrainingSet. It is not his actual training set, which I did not find in print, just my simulation.

Make a plot like Bishop’s figure 1.7 (page 10).

Partials: Gradients of the Unknown Parameters

The Observation Equations

MLE: The Normal Equations

For comparison purposes, mechanize the normal equations. Expect them to over-fit. Here is an example for order 3.

The normal equations as a symbolic polynomial on bts follows. Notice we can increase the order beyond the number of data, creating an underdetermined system. Underdetermined systems are not typical in real-world data processing. Usually the number of data exceed the order and the system is overdetermined. The pseudoinverse in mleFit is agnostic to the distinction.

RLS: Recurrent Least Squares

Non-Renormalized RLS Fitting Procedure

Here is an example of a polynomial of order three: four parameters including bias.

KAL: Foldable Kalman Filter

The foldable Kalman filter (KAL) follows below. This version has only the update phase of a typical Kalman filter because the parameters-to-estimate are constant and therefore there is no predict phase.

Kalman Fitting Procedure

See All Three

The following interactive demonstration shows mleFit (normal equations), rlsFit (non-renormalized recurrent least squares), and kalFit (Kalman folding) on Bishop’s training set. You may choose the order of the polynomial to fit and the logarithms of the a-priori information matrix (for RLE) and the covariance (for KAL).

Run the polynomial order up to nine, then run -logΛ and logσξ2 all the way to the right, to their maximum values, to see overfitting, that is, failure of regularization.

The Green line is Ground Truth.

◼

Renormalizing RLS

The following slogan captures regularization in a nutshell.

As observation noise decreases, the solutions believe the observations more than the a-prioris and RLS and KAL over-fit. As a-priori covariance decreases, RLS and KAL believe the a-prioris more than the observations and the solution regularizes.

Renormalized RLS

RLS, so far, is normalized to unit observation noise. How to modify RLS to account for non-unit observation noise?

Redefine rlsUpdate to include this renormalization and notice that it’s always equal to KAL.

Renormalized RLS beats KAL

KAL and renormalized RLS are mathematically equivalent (stated without proof). Operationally, KAL subtracts matrices to update estimated covariance. Subtraction is exposed to catastrophic cancelation. Renormalized RLS only adds to the information matrix, so is exposed only to ill-conditioning, which is empirically less severe than catastrophic cancelation. We show this below.

Regularization and MAP

◼

Bishop’s MAP

The MAP Equations

Here are some links between Bishop’s formulation and ours, without derivation.

Φ Vectors

S Inverse

Bishop’s Equation 1.72.

MAP Mean

Bishop’s Equation 1.70, replacing Inverse with LinearSolve as we usually do:

a-priori Variances α and β

Mean Is Invariant Under “Swap and Invert”

Numerical Equivalences

RLS does not suffer from these ills because it never subtracts. RLS is still exposed to ill-conditioning of the information matrix, but that seems numerically less harmful to the final result in this example.

Wrap RLS in Quiet to suppress warnings. There is no free lunch; MAP also shows ill-conditioning and is similarly wrapped.

Recurrent RLS Fitting Procedure

Covariance and Information Matrices

Covariance of the Prediction

Policy Gradient is co-RLS, co-KAL

straight transposes of Equation 3. LinearSolve operates on the transposed right-hand side of the recurrence. Transpose the solution to get the recurrence. Apply the resulting dual model to the transpose of the original data. The following function is the dual of RLS update, and the fold effects co-RLS. We check it against our running slope-intercept example. We have shown that RLS is equivalent to KAL, thus co-RLS is equivalent to co-KAL.

◼

Finite-Difference Policy Gradient

The finite-difference method of policy-gradient is an instance of co-RLS (see http://www.scholarpedia.org/article/Policy_gradient_methods).

The Moore-Penrose right pseudoinverse RPI

furnishes a large-memory, slow, numerically hazardous solution of the problem:

Padé Approximant via NIST

A Padé approximant is a ratio of two polynomials, where the bias term in the denominator is unity to prevent division by zero. Quoting from Srini Kumar and Bob Horton https://blog.revolutionanalytics.com/2017/04/fitting-rational-functions-with-lm.html:

The following example from NIST (https://www.itl.nist.gov/div898/strd/nls/data/LINKS/DATA/Thurber.dat, edited to remove some blank lines) lets us illustrate:

NIST/ITL StRD
Dataset Name: Thurber (Thurber.dat)

File Format: ASCII
Starting Values (lines 41 to 47)
Certified Values (lines 41 to 52)
Data (lines 61 to 97)

Procedure: Nonlinear Least Squares Regression

Description: These data are the result of a NIST study involving
semiconductor electron mobility. The response
variable is a measure of electron mobility, and the
predictor variable is the natural log of the density.

Reference: Thurber, R., NIST (197?).
Semiconductor electron mobility modeling.

Data: 1 Response Variable (y = electron mobility)
1 Predictor Variable (x = log[density])
37 Observations
Higher Level of Difficulty
Observed Data

Model: Rational Class (cubic/cubic)
7 Parameters (b1 to b7)

y = (b1 + b2*x + b3*x**2 + b4*x**3) /
(1 + b5*x + b6*x**2 + b7*x**3) + e

Starting Values Certified Values

Start 1 Start 2 Parameter Standard Deviation
b1 = 1000 1300 1.2881396800E+03 4.6647963344E+00
b2 = 1000 1500 1.4910792535E+03 3.9571156086E+01
b3 = 400 500 5.8323836877E+02 2.8698696102E+01
b4 = 40 75 7.5416644291E+01 5.5675370270E+00
b5 = 0.7 1 9.6629502864E-01 3.1333340687E-02
b6 = 0.3 0.4 3.9797285797E-01 1.4984928198E-02
b7 = 0.03 0.05 4.9727297349E-02 6.5842344623E-03

Residual Sum of Squares: 5.6427082397E+03
Residual Standard Deviation: 1.3714600784E+01
Degrees of Freedom: 30
Number of Observations: 37

Data: y x
80.574E0 -3.067E0
84.248E0 -2.981E0
87.264E0 -2.921E0
87.195E0 -2.912E0
89.076E0 -2.840E0
89.608E0 -2.797E0
89.868E0 -2.702E0
90.101E0 -2.699E0
92.405E0 -2.633E0
95.854E0 -2.481E0
100.696E0 -2.363E0
101.060E0 -2.322E0
401.672E0 -1.501E0
390.724E0 -1.460E0
567.534E0 -1.274E0
635.316E0 -1.212E0
733.054E0 -1.100E0
759.087E0 -1.046E0
894.206E0 -0.915E0
990.785E0 -0.714E0
1090.109E0 -0.566E0
1080.914E0 -0.545E0
1122.643E0 -0.400E0
1178.351E0 -0.309E0
1260.531E0 -0.109E0
1273.514E0 -0.103E0
1288.339E0 0.010E0
1327.543E0 0.119E0
1353.863E0 0.377E0
1414.509E0 0.790E0
1425.208E0 0.963E0
1421.384E0 1.006E0
1442.962E0 1.115E0
1464.350E0 1.572E0
1468.705E0 1.841E0
1447.894E0 2.047E0
1457.628E0 2.200E0

Read the Data into Mathematica

Symbolic Model

With the notation of the NIST example, rather than that of Kumar and Horton:

Linearized Model

Residual

Partials

Priors

Define global ξ01 as “Start1” priors in the NIST document:

Define global ξ02 as “Start2” in the NIST document:

"Certified Values" in the NIST document :

Input collection for recurrent form:

Setup for RLS

Define a symbolic column vector of coefficients to estimate. Set up some arbitrary prior covariances. Set up rules for evaluating symbolic coefficients at given numerical values for plotting.

TODO: Fix RLS so it can handle this scenario.

As usual, the observation noise, sqrtPΖ, is a constant over the run.

Unit test of rlsUpdate:

It turns out that the RLS update is slow and uninstructive. kalmanUpdate suffices to solve the problem.

We leave detailed examination of the covariances of this solution to future work.

Conclusion

Kalman folding (KAL) produces the same results as renormalized recurrent least squares (RLS) and as maximum a-posteriori (MAP) for appropriate choices of KAL and RLS covariances, or regularization hyperparameters for MAP.

Numerically and symbolically through order 3, MAP produces the same estimates as KAL and RLS, though not covariances, when its hyperparameters are swapped and inverted. MAP covariances are recovered in post-processing by un-swapping and re-inverting.

Finite-differences for policy gradient (FDPG) in reinforcement learning is co-RLS, co-KAL.

KAL and RLS have much better space-time efficiency than non-recurrent MAP. They avoid storing and multiplying large matrices. In all cases, we avoid matrix inverses by solving linear systems internally.

CITE THIS NOTEBOOK

Efficient Bayesian Regularization
by Brian Beckman
Wolfram Community, STAFF PICKS, January 22, 2024
https://community.wolfram.com/groups/-/m/t/3107102