<figure><img src="image_1.png" /></figure>

HAL
open science

# Revisiting Whittaker-Henderson Smoothing

Guillaume Biessy

► To cite this version:

Guillaume Biessy. Revisiting Whittaker-Henderson Smoothing. 2023. hal-04124043v2

HAL Id: hal-04124043

https://hal.science/hal-04124043v2

Preprint submitted on 18 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Whittaker–Henderson smoothing revisited: A modern statistical framework for practical use

Guillaume Biessy

► To cite this version:

Guillaume Biessy. Whittaker–Henderson smoothing revisited: A modern statistical framework for practical use. ASTIN Bulletin, In press, 10.1017/asb.2025.10061. hal-04124043v5

HAL Id: hal-04124043

https://hal.science/hal-04124043v5

Submitted on 3 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Revisiting Whittaker-Henderson Smoothing

Guillaume Biessy*, PhD, LinkPact†and Sorbonne Université‡

September 16, 2023

Introduced nearly a century ago, Whittaker-Henderson smoothing remains one of the most commonly used methods by actuaries for constructing one-dimensional and two-dimensional experience tables for mortality and other Life Insurance risks. This paper proposes to reframe this smoothing technique within a modern statistical framework and addresses six questions of practical interest regarding its use.

Firstly, we adopt a Bayesian view of this smoothing method to build credible intervals. Next, we shed light on the choice of observation vectors and weights to which the smoothing should be applied by linking it to a maximum likelihood estimator introduced in the context of duration models. We then enhance the precision of the smoothing by relaxing an implicit asymptotic approximation on which it relies. Afterward, we select the smoothing parameters based on maximizing a marginal likelihood. We later improve numerical performance in the presence of a large number of observation points and, consequently, parameters. Finally, we extrapolate the results of the smoothing while preserving consistency between estimated and predicted values through the use of constraints.

**Keywords:** smoothing methods, duration models, experience tables, maximum likelihood, generalized additive models, empirical Bayesian approach, marginal likelihood, extrapolation.

# 1 Introduction

## 1.1 Notations

In this paper, vector names are written in bold characters and matrix names in uppercase letters. If $\mathbf{y}$ is a vector and $\mathbf{A}$ is a matrix, $\mathrm{Var}(\mathbf{y})$ denotes the variance-covariance matrix

<footnote>*guillaume.biessy78@gmail.com</footnote>

<footnote>†LinkPact, 80 boulevard Pasteur, 75015 Paris, France</footnote>

<footnote>‡Sorbonne Université, CNRS, Laboratoire de Probabilités, Statistique et Modélisation, LPSM, 4 place Jussieu, 75005 Paris, France</footnote>
# Whittaker-Henderson Smoothing Revisited: A Modern Statistical Framework for Practical Use

Guillaume Biessy*, PhD, LinkPact†and Sorbonne Université‡

September 3, 2025

Introduced over a century ago, Whittaker-Henderson smoothing remains widely used by actuaries in constructing one-dimensional and two-dimensional experience tables for mortality, disability and other life insurance risks. In this paper, we reinterpret this smoothing technique within a modern statistical framework and address six practically relevant questions about its use.

First, we adopt a Bayesian perspective on this method to construct credible intervals. Second, in the context of survival analysis, we clarify how to choose the observation and weight vectors by linking the smoothing technique to a maximum likelihood estimator. Third, we improve accuracy by relaxing the method's reliance on an implicit normal approximation. Fourth, we select the smoothing parameters by maximizing a marginal likelihood function. Fifth, we improve computational efficiency when dealing with numerous observation points and consequently parameters. Finally, we develop an extrapolation procedure that ensures consistency between estimated and predicted values through constraints.

## Table of contents

<table><tr><td><b>1 Introduction</b></td><td><b>3</b></td></tr><tr><td>1.1 A brief reminder of WH smoothing mathematical formulation</td><td>3</td></tr><tr><td>1.2 Structure of the paper</td><td>5</td></tr><tr><td><b>2 How to measure uncertainty in smoothing results?</b></td><td><b>9</b></td></tr><tr><td>2.1 Maximum a posteriori estimate</td><td>9</td></tr></table>

<footnote>*guillaume.biessy78@gmail.com</footnote>

<footnote>†LinkPact, 75015 Paris, France</footnote>

<footnote>‡Sorbonne Université, CNRS, Laboratoire de Probabilités, Statistique et Modélisation, LPSM, 75005 Paris, France</footnote>
associated with **y**, **diag**(A) represents the diagonal of matrix A, and Diag(**y**) is the diagonal matrix such that **diag**(Diag(**y**)) = **y**. The sum of the diagonal elements of A is denoted as tr(A) and its transpose as A^T. In the case where A is invertible, A^(-1) denotes its inverse, A^(-T) is the inverse of its transpose, and |A| is the product of the eigenvalues of A. For a non-invertible matrix A, A^(-) refers to the Moore-Penrose pseudoinverse of A, and |A|_+ denotes the product of the non-zero eigenvalues of A. By writing the eigendecomposition as A = UΣV^T, where U and V are two orthogonal matrices and Σ is a diagonal matrix containing the eigenvalues of A, and by denoting Σ^(-) as the matrix obtained by replacing the non-zero eigenvalues in Σ with their inverses while keeping the zero eigenvalues unchanged, the pseudoinverse is given by A^(-) = VΣ^(-)U^T. The Kronecker product of two matrices A and B is denoted as A ⊗ B, and their Hadamard product, or element-wise product, is denoted as A ⊙ B. **vec**(A) represents the vector obtained by stacking the columns of matrix A together. Finally, the symbol ∝ signifies a proportional relationship between the expressions on both sides of it.

## 1.2 Origin

Whittaker-Henderson smoothing is a gradation method aimed at correcting the effect of sampling fluctuations on an observation vector. It is applied to evenly-spaced discrete observations. Initially proposed by Whittaker (1922) for constructing mortality tables and further developed by the works of Henderson (1924), it remains one of the most popular methods among actuaries for constructing experience tables in life insurance. Extending to two-dimensional tables, it can be used for studying various risks, including but not limited to mortality, disability, long-term care, lapse, mortgage default, and unemployment.

## 1.3 The one-dimensional case

Let **y** be a vector of observations and **w** a vector of positive weights, both of size *n*. The estimator associated with Whittaker-Henderson smoothing is given by:

$$ \hat{\mathbf{y}} = \underset{\theta}{\operatorname{argmin}} \{F(\mathbf{y}, \mathbf{w}, \theta) + R_{\lambda,q}(\theta)\} \quad (1) $$

where:

* $F(\mathbf{y}, \mathbf{w}, \theta) = \sum_{i=1}^{n} w_i (y_i - \theta_i)^2$ represents a fidelity criterion to the observations,

* $R_{\lambda,q}(\theta) = \lambda \sum_{i=1}^{n-q} (\Delta^q \theta)_i^2$ represents a smoothness criterion.
<table>
  <tbody>
    <tr>
      <td>2.2</td>
      <td>Posterior distribution of θ | y</td>
      <td>9</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Consequence for the WH smoothing</td>
      <td>10</td>
    </tr>
    <tr>
      <td><b>3</b></td>
      <td><b>Which observation and weight vectors to use?</b></td>
      <td><b>10</b></td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Survival analysis framework</td>
      <td>10</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Likelihood equations</td>
      <td>11</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Consequence for the WH smoothing</td>
      <td>12</td>
    </tr>
    <tr>
      <td><b>4</b></td>
      <td><b>How to improve the accuracy of smoothing with limited data volume?</b></td>
      <td><b>12</b></td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Generalized Whittaker-Henderson smoothing</td>
      <td>12</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Impact of the normal approximation in the original smoothing</td>
      <td>13</td>
    </tr>
    <tr>
      <td><b>5</b></td>
      <td><b>How to select the smoothing parameters?</b></td>
      <td><b>15</b></td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Impact of smoothing parameter choice</td>
      <td>15</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Statistical criteria for parameter selection</td>
      <td>16</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Selection in the original smoothing</td>
      <td>18</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Selection in the generalized smoothing</td>
      <td>19</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Algorithms for the maximization of the marginal likelihood</td>
      <td>20</td>
    </tr>
    <tr>
      <td>5.6</td>
      <td>Performance comparison</td>
      <td>21</td>
    </tr>
    <tr>
      <td><b>6</b></td>
      <td><b>How to improve smoothing computational efficiency?</b></td>
      <td><b>23</b></td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Motivation</td>
      <td>23</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Practical computation for penalized smoothers</td>
      <td>25</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Banded optimization for WH smoothing</td>
      <td>27</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Natural parameterization and rank reduction of WH smoothing</td>
      <td>29</td>
    </tr>
    <tr>
      <td><b>7</b></td>
      <td><b>How to extrapolate the smoothing?</b></td>
      <td><b>34</b></td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Defining the extrapolation of the smoothing</td>
      <td>34</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Unconstrained solution for the 1D case</td>
      <td>36</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Constrained solution for the 2D case</td>
      <td>37</td>
    </tr>
    <tr>
      <td><b>8</b></td>
      <td><b>Discussion</b></td>
      <td><b>40</b></td>
    </tr>
    <tr>
      <td><b>References</b></td>
      <td></td>
      <td><b>43</b></td>
    </tr>
    <tr>
      <td><b>Appendices</b></td>
      <td></td>
      <td><b>46</b></td>
    </tr>
    <tr>
      <td>A</td>
      <td>Exposure computation in the survival analysis framework</td>
      <td>46</td>
    </tr>
    <tr>
      <td>B</td>
      <td>Simulated datasets</td>
      <td>48</td>
    </tr>
    <tr>
      <td>C</td>
      <td>Derivation of marginal likelihood and LAML</td>
      <td>52</td>
    </tr>
    <tr>
      <td>D</td>
      <td>Algorithms</td>
      <td>54</td>
    </tr>
    <tr>
      <td>E</td>
      <td>Illustrations for the natural parameterization of WH smoothing</td>
      <td>56</td>
    </tr>
    <tr>
      <td>F</td>
      <td>Derivation of constrained extrapolation in the 2D case</td>
      <td>60</td>
    </tr>
    <tr>
      <td>G</td>
      <td>Choosing the order of the difference matrices</td>
      <td>62</td>
    </tr>
  </tbody>
</table>
In the latter expression, $\Delta^q$ denotes the forward difference operator of order $q$, such that for any $i \in [1, n-q]$:

$$ (\Delta^q \boldsymbol{\theta})_i = \sum_{k=0}^{q} \binom{q}{k} (-1)^{q-k} \theta_{i+k}. $$

Let us define $W = \text{Diag}(\mathbf{w})$, the diagonal matrix of weights, and $D_{n,q}$ as the order $q$ difference matrix of dimensions $(n-q) \times n$, such that $(D_{n,q}\boldsymbol{\theta})_i = (\Delta^q\boldsymbol{\theta})_i$ for all $i \in [1, n-q]$. The most commonly used difference matrices of order 1 and 2 have the following forms:

$$ D_{n,1} = \begin{bmatrix} -1 & 1 & 0 & \dots & 0 \\ 0 & -1 & 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & -1 & 1 \end{bmatrix} \quad \text{and} \quad D_{n,2} = \begin{bmatrix} 1 & -2 & 1 & 0 & \dots & 0 \\ 0 & 1 & -2 & 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & \dots & 0 & 1 & -2 & 1 \end{bmatrix}. $$

The fidelity and smoothness criteria can be rewritten in matrix form as:

$$ F(\mathbf{y}, \mathbf{w}, \boldsymbol{\theta}) = (\mathbf{y} - \boldsymbol{\theta})^T W (\mathbf{y} - \boldsymbol{\theta}) \quad \text{and} \quad R_{\lambda,q}(\boldsymbol{\theta}) = \lambda \boldsymbol{\theta}^T D_{n,q}^T D_{n,q} \boldsymbol{\theta}. $$

The associated estimator for smoothing becomes:

$$ \hat{\mathbf{y}} = \underset{\boldsymbol{\theta}}{\operatorname{argmin}} \left\{ (\mathbf{y} - \boldsymbol{\theta})^T W (\mathbf{y} - \boldsymbol{\theta}) + \boldsymbol{\theta}^T P_{\lambda} \boldsymbol{\theta} \right\} \quad (2) $$

where $P_{\lambda} = \lambda D_{n,q}^T D_{n,q}$.

## 1.4 The two-dimensional case

In the two-dimensional case, let us consider a matrix $Y$ of observations and a matrix $\Omega$ of non-negative weights, both of dimensions $n_x \times n_z$. The estimator associated with the Whittaker-Henderson smoothing can be written as:

$$ \hat{Y} = \underset{\Theta}{\operatorname{argmin}} \{F(Y, \Omega, \Theta) + R_{\lambda,q}(\Theta)\} $$

where:

* $F(Y, \Omega, \Theta) = \sum_{i=1}^{n_x} \sum_{j=1}^{n_z} \Omega_{i,j} (Y_{i,j} - \Theta_{i,j})^2$ represents a fidelity criterion to the observations,

* $R_{\lambda,q}(\Theta) = \lambda_x \sum_{j=1}^{n_z} \sum_{i=1}^{n_x} (\Delta^{q_x} \Theta_{\bullet,j})_i^2 + \lambda_z \sum_{i=1}^{n_x} \sum_{j=1}^{n_z} (\Delta^{q_z} \Theta_{i,\bullet})_j^2$ is a smoothness criterion.
# Notations

In this paper, vectors are denoted in boldface and matrix names in uppercase letters. If $\mathbf{y}$ is a vector and $A$ is a matrix, $\mathrm{Var}(\mathbf{y})$ denotes the variance-covariance matrix associated with $\mathbf{y}$, $\mathrm{diag}(A)$ represents the diagonal of matrix $A$, and $\mathrm{Diag}(\mathbf{y})$ is the diagonal matrix such that $\mathrm{diag}(\mathrm{Diag}(\mathbf{y})) = \mathbf{y}$. The sum of the diagonal elements of $A$ is denoted as $\mathrm{tr}(A)$ and its transpose as $A^T$. In the case where $A$ is invertible, $A^{-1}$ denotes its inverse and $|A|$ denotes the product of the eigenvalues of $A$. For a non-invertible matrix $A$, $A^{-1}$ refers to the Moore-Penrose pseudo-inverse of $A$, and $|A|_+$ denotes the product of the non-zero eigenvalues of $A$. By writing the eigendecomposition as $A = U\Sigma V^T$, where $U$ and $V$ are orthogonal matrices and $\Sigma$ is a diagonal matrix containing the eigenvalues of $A$, and by denoting $\Sigma^{-1}$ as the matrix obtained by replacing the non-zero eigenvalues in $\Sigma$ with their inverses leaving the zero eigenvalues unchanged, the pseudo-inverse is given by $A^{-1} = V\Sigma^{-1}U^T$. The Kronecker product of two matrices $A$ and $B$ is denoted as $A \otimes B$, and their Hadamard (element-wise) product, is denoted as $A \odot B$. $\lfloor x \rfloor$ denotes the greatest integer less than or equal to $x \in \mathbb{R}$. Finally, the symbol $\propto$ denotes proportionality between the expressions on both sides.

# 1 Introduction

Whittaker-Henderson (WH) smoothing is a graduation method designed to mitigate the effects of sampling fluctuations in a vector of evenly spaced discrete observations. Although this method was originally proposed by Bohlmann (1899), it is named after Whittaker (1923), who applied it to graduate mortality tables, and Henderson (1924), who popularized it among actuaries in the United States. The method was later extended to two dimensions by Knorr (1984). WH smoothing may be used to build experience tables for a broad spectrum of life insurance risks, such as mortality, disability, long-term care, lapse, mortgage default and unemployment. We begin with a brief overview of the method before outlining the structure and main contributions of the paper.

## 1.1 A brief reminder of WH smoothing mathematical formulation

### The one-dimensional case

Let $\mathbf{y}$ be a vector of observations and $\mathbf{w}$ a vector of positive weights, both of size $n$. The estimator associated with Whittaker-Henderson smoothing is given by:

$$ \hat{\mathbf{y}} = \underset{\theta}{\operatorname{argmin}} \{F(\mathbf{y}, \mathbf{w}, \theta) + R_{\lambda,q}(\theta)\} \quad (1) $$

where:
This latter criterion can be written as the sum of two one-dimensional regularization criteria, with orders $q_x$ and $q_z$, applied respectively to all rows and all columns of $\Theta$. It is also possible to adopt matrix notations by defining $\mathbf{y} = \text{vec}(Y)$, $\mathbf{w} = \text{vec}(\Omega)$, and $\boldsymbol{\theta} = \text{vec}(\Theta)$ as the vectors obtained by stacking the columns of the matrices $Y$, $\Omega$, and $\Theta$, respectively. Additionally, let us denote $W = \text{Diag}(\mathbf{w})$ and $n = n_x \times n_z$. The fidelity and smoothness criteria can then be rewritten as linear combinations of the vectors $\mathbf{y}$, $\mathbf{w}$, and $\boldsymbol{\theta}$:

$$F(\mathbf{y}, \mathbf{w}, \boldsymbol{\theta}) = (\mathbf{y} - \boldsymbol{\theta})^T W (\mathbf{y} - \boldsymbol{\theta}) \\ R_{\lambda,q}(\boldsymbol{\theta}) = \boldsymbol{\theta}^T (\lambda_x I_{n_z} \otimes D_{n_x,q_x}^T D_{n_x,q_x} + \lambda_z D_{n_z,q_z}^T D_{n_z,q_z} \otimes I_{n_x}) \boldsymbol{\theta}.$$

and the estimator associated with the smoothing takes the form:

$$\hat{\mathbf{y}} = \underset{\boldsymbol{\theta}}{\operatorname{argmin}} \left\{ (\mathbf{y} - \boldsymbol{\theta})^T W (\mathbf{y} - \boldsymbol{\theta}) + \boldsymbol{\theta}^T P_{\lambda} \boldsymbol{\theta} \right\}$$

of Equation 2 except in this case $P_{\lambda} = \lambda_x I_{n_z} \otimes D_{n_x,q_x}^T D_{n_x,q_x} + \lambda_z D_{n_z,q_z}^T D_{n_z,q_z} \otimes I_{n_x}$.

## 1.5 Explicit solution

Equation 2 has an explicit solution:

$$\hat{\mathbf{y}} = (W + P_{\lambda})^{-1} W \mathbf{y}. \qquad (3)$$

Indeed, as a minimum, $\hat{\mathbf{y}}$ satisfies:

$$0 = \frac{\partial}{\partial \boldsymbol{\theta}} \Big|_{\hat{\mathbf{y}}} \left\{ (\mathbf{y} - \boldsymbol{\theta})^T W (\mathbf{y} - \boldsymbol{\theta}) + \boldsymbol{\theta}^T P_{\lambda} \boldsymbol{\theta} \right\} = -2X^T W (\mathbf{y} - \hat{\mathbf{y}}) + 2X^T P_{\lambda} \hat{\mathbf{y}}.$$

It follows that $(W + P_{\lambda})\hat{\mathbf{y}} = W\mathbf{y}$, and if $W + P_{\lambda}$ is invertible, then $\hat{\mathbf{y}}$ is indeed a solution of Equation 3. When $\lambda \neq 0$, $W + P_{\lambda}$ is invertible as long as $\mathbf{w}$ has $q$ non-zero elements in the one-dimensional case, and $\Omega$ has at least $q_x \times q_z$ non-zero elements distributed over $q_x$ different rows and $q_z$ different columns in the two-dimensional case. These sufficient conditions are always satisfied in practice and will not be demonstrated here.
• $F(\mathbf{y}, \mathbf{w}, \boldsymbol{\theta}) = \sum_{i=1}^{n} w_i (y_i - \theta_i)^2$ represents a fidelity criterion with respect to the observations,

• $R_{\lambda,q}(\boldsymbol{\theta}) = \lambda \sum_{i=1}^{n-q} (\Delta^q \boldsymbol{\theta})_i^2$ represents a smoothness criterion.

In the latter expression, $\lambda \ge 0$ is a smoothing parameter and $\Delta^q$ denotes the forward difference operator of order $q$, such that for any $i \in \{1, \dots, n-q\}$:

$$ (\Delta^q \boldsymbol{\theta})_i = \sum_{k=0}^{q} \binom{q}{k} (-1)^{q-k} \theta_{i+k}. $$

Define $W = \text{Diag}(\mathbf{w})$, the diagonal matrix of weights, and $D_{n,q}$ as the order $q$ difference matrix of dimensions $(n-q) \times n$, such that $(D_{n,q}\boldsymbol{\theta})_i = (\Delta^q \boldsymbol{\theta})_i$ for all $i \in [1, n-q]$. The first- and second-order difference matrices are given by:

$$ D_{n,1} = \begin{bmatrix} -1 & 1 & 0 & \dots & 0 \\ 0 & -1 & 1 & \dots & \vdots \\ \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \dots & 0 & -1 & 1 \end{bmatrix} \quad \text{and} \quad D_{n,2} = \begin{bmatrix} 1 & -2 & 1 & 0 & \dots & 0 \\ 0 & 1 & -2 & 1 & \dots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & \dots & 0 & 1 & -2 & 1 \end{bmatrix}. $$

while higher-order difference matrices follow the recursive formula $D_{n,q} = D_{n-1,q-1}D_{n,1}$. The fidelity and smoothness criteria can be rewritten with matrix notations as:

$$ F(\mathbf{y}, \mathbf{w}, \boldsymbol{\theta}) = (\mathbf{y} - \boldsymbol{\theta})^T W (\mathbf{y} - \boldsymbol{\theta}) \quad \text{and} \quad R_{\lambda,q}(\boldsymbol{\theta}) = \lambda \boldsymbol{\theta}^T D_{n,q}^T D_{n,q} \boldsymbol{\theta} $$

and the WH smoothing estimator thus becomes:

$$ \hat{\mathbf{y}} = \underset{\boldsymbol{\theta}}{\operatorname{argmin}} \left\{ (\mathbf{y} - \boldsymbol{\theta})^T W (\mathbf{y} - \boldsymbol{\theta}) + \boldsymbol{\theta}^T P_{\lambda} \boldsymbol{\theta} \right\} \quad (2) $$

where $P_{\lambda} = \lambda D_{n,q}^T D_{n,q}$.

## The two-dimensional case

In the two-dimensional case, consider a matrix $Y$ of observations and a matrix $\Omega$ of non-negative weights, both of dimensions $n_x \times n_z$. The WH smoothing estimator solves:

$$ \hat{Y} = \underset{\Theta}{\operatorname{argmin}} \{F(Y, \Omega, \Theta) + R_{\lambda,q}(\Theta)\} $$

where:

• $F(Y, \Omega, \Theta) = \sum_{i=1}^{n_x} \sum_{j=1}^{n_z} \Omega_{i,j} (Y_{i,j} - \Theta_{i,j})^2$ represents a fidelity criterion with respect to the observations,
1.6 Impact of the smoothing parameter(s)

In Equation 1 (y – θ)<sup>T</sup>W(y – θ) represents a fidelity criterion to the observations and θ<sup>T</sup>P<sub>λ</sub>θ represents a smoothness criterion. The relative importance of these criteria is controlled by the parameter (or pair of parameters in the two-dimensional case) λ. Figure 1 provides an illustration of one-dimensional smoothing on mortality data using three different values of the smoothing parameter. The effective degrees of freedom shown in this figure can be calculated by summing the diagonal values of H, the *hat matrix* of the model. The hat matrix satisfies ŷ = Hy and can be identified in Equation 3 as H = (W + P<sub>λ</sub>)<sup>-1</sup>W. These effective degrees of freedom serve as a non-parametric equivalent of the number of independent parameters in parametric models but can take non-integer values. The concept of degrees of freedom will be discussed in more detail in Section 6 of the paper. It can be observed here that the result of smoothing is highly sensitive to the chosen value of the smoothing parameter. The choice of λ = 10<sup>1</sup> leads to a highly volatile fit that largely reproduces the sampling fluctuations present in the data. On the other hand, the choice of λ = 10<sup>7</sup> seems too rigid to capture the underlying pattern of the data. The choice of λ = 10<sup>4</sup> appears, at first glance, as a satisfactory compromise between these two extremes.

<figure><img src="image_3.png" /><figcaption>Figure 1: Whittaker-Henderson smoothing applied to a portfolio of synthetic mortality data for 3 choices of the smoothing parameter. The points represent the initial observations, the curves depict the result of the smoothing, and the shaded areas represent the associated credibility intervals. *edf* refers to the effective degrees of freedom of the model.</figcaption></figure>

In the two-dimensional case, two penalties are applied to the model parameters. Figure 2 illustrates the application of smoothing to mortality data for a portfolio of insured individuals with dependence, based on the age of the insured individuals and their duration in the dependent state. Different combinations of the two smoothing parameters in this model are used to produce the 9 smoothers represented in this figure. While some combinations of smoothing parameters
• $R_{\lambda,q}(\Theta) = \lambda_x \sum_{j=1}^{n_z} \sum_{i=1}^{n_x-q_x} (\Delta^{q_x}\Theta_{\bullet,j})_i^2 + \lambda_z \sum_{i=1}^{n_x} \sum_{j=1}^{n_z-q_z} (\Delta^{q_z}\Theta_{i,\bullet})_j^2$ is a smoothness criterion with $\lambda = (\lambda_x, \lambda_z)$.

This latter criterion adds row-wise and column-wise regularization criteria to $\Theta$, with respective orders $q_x$ and $q_z$, weighted by non-negative smoothing parameters $\lambda_x$ and $\lambda_z$. In matrix notation, let $\mathbf{y} = \text{vec}(\mathbf{Y})$, $\mathbf{w} = \text{vec}(\Omega)$, and $\boldsymbol{\theta} = \text{vec}(\Theta)$ as the vectors obtained by stacking the columns of the matrices $\mathbf{Y}$, $\Omega$, and $\Theta$, respectively. Additionally, denote $\mathbf{W} = \text{Diag}(\mathbf{w})$ and $n = n_x \times n_z$. The fidelity and smoothness criteria become:

$$F(\mathbf{y}, \mathbf{w}, \boldsymbol{\theta}) = (\mathbf{y} - \boldsymbol{\theta})^T \mathbf{W} (\mathbf{y} - \boldsymbol{\theta}) \\ R_{\lambda,q}(\boldsymbol{\theta}) = \boldsymbol{\theta}^T (\lambda_x I_{n_z} \otimes D_{n_x,q_x}^T D_{n_x,q_x} + \lambda_z D_{n_z,q_z}^T D_{n_z,q_z} \otimes I_{n_x}) \boldsymbol{\theta}$$

and the associated estimator also takes the form of Equation 2 except in this case

$$P_{\lambda} = \lambda_{x} I_{n_{z}} \otimes D_{n_{x},q_{x}}^{T} D_{n_{x},q_{x}} + \lambda_{z} D_{n_{z},q_{z}}^{T} D_{n_{z},q_{z}} \otimes I_{n_{x}}.$$

Extension to higher dimensions is straightforward and not discussed here.

## An explicit solution

If $W + P_\lambda$ is invertible, Equation 2 admits the closed-form solution:

$$\hat{\mathbf{y}} = (W + P_{\lambda})^{-1} W \mathbf{y}. \qquad (3)$$

Indeed, as a minimum, $\hat{\mathbf{y}}$ satisfies:

$$0 = \left. \frac{\partial}{\partial \boldsymbol{\theta}} \right|_{\hat{\mathbf{y}}} \{ (\mathbf{y} - \boldsymbol{\theta})^T W (\mathbf{y} - \boldsymbol{\theta}) + \boldsymbol{\theta}^T P_{\lambda} \boldsymbol{\theta} \} = -2W(y - \hat{\mathbf{y}}) + 2P_{\lambda} \hat{\mathbf{y}}.$$

It follows that $(W + P_\lambda)\hat{\mathbf{y}} = W\mathbf{y}$, proving Equation 3. If $\lambda \neq 0$, $W + P_\lambda$ is invertible as long as $\mathbf{w}$ has $q$ non-zero elements in the one-dimensional case, and $\Omega$ has at least $q_x \times q_z$ non-zero elements spread across $q_x$ different rows and $q_z$ different columns in the two-dimensional case. These conditions are always met in real datasets.

## 1.2 Structure of the paper

Introduced a century ago, Whittaker-Henderson (WH) smoothing remains widely used by actuaries, particularly in France and North America (Canadian Institute of Actuaries 2017; Society of Actuaries 2018). Other non-parametric smoothing methods have since emerged, notably spline-based techniques (Reinsch 1967), which gained even greater popularity with P-splines (Eilers and Marx 1996). A broader overview of alternative smoothers is available in Wood (2017, chap. 5).
may be deemed implausible, such as those associated with the choice $\lambda_x = 10^0$, it is challenging in the two-dimensional case to determine graphically which combination of parameters is most suitable. The selection of the smoothing parameter will be discussed in detail in Section 5.

<figure><img src="image_4.png" /><figcaption>Figure 2: Whittaker-Henderson smoothing applied to a synthetic portfolio of policyholders in long-term care situation for 9 combinations of smoothing parameters. The contour lines and colors depict the surface of the smoothed mortality force as a function of the age of the policyholders and the duration in the long-term care state.</figcaption></figure>

## 1.7 Plan for the paper

This paper aims to address 6 practical questions associated with the use of Whittaker-Henderson smoothing, each of which will be covered in a dedicated section.
For evenly spaced discrete observations, WH smoothing may be considered a particular case of P-splines with degree-zero splines and identity model matrix. Its appeal lies in its simplicity: no selection of knots, parameters equal to fitted values, and shape controlled solely via penalization. However, it involves more parameters than low-rank smoothers, making it more computationally intensive.

Originally proposed as an empirical alternative to polynomial regression and weighted averages, WH smoothing offered key benefits noted by Whittaker (1923): first *q* moment preservation, adjustable smoothing parameters, and robustness at boundaries. While smoothing theory has evolved—particularly via generalized additive models (Hastie and Tibshirani 1990), use of WH smoothing by actuaries remains largely unchanged. This paper reinterprets WH within modern statistical theory to bridge that gap and address six practical questions, each discussed in a dedicated section.

**How to measure uncertainty in smoothing results?**

We propose a method to quantify the uncertainty in WH smoothing based on data volume, a topic that has received little attention in the literature. In a Frequentist framework, the WH estimator is biased, which complicates the construction of valid confidence intervals for finite samples. However, under certain conditions, WH smoothing can be viewed as a Bayesian model, enabling the derivation of credible intervals. This Bayesian interpretation was originally suggested by Whittaker (1923) as a justification for the method and formally revisited decades later by Taylor (1992). In this section, we build on that equivalence to derive credible intervals for WH smoothing.

**Which observation and weight vectors to use?**

For the Bayesian interpretation of WH smoothing discussed in Section 2 to hold, it must be applied to a vector **y** of independent, normally distributed observations with known variances. The weight vector **w** should then contain the inverse variances (up to a constant), as noted by Taylor (1992) and Verrall (1993). We show that, under piecewise constant transition intensities in duration models, the maximum likelihood estimator of crude rates produces vectors (**y**, **w**) that asymptotically meet these conditions. This, combined with the results from the previous section, offers a statistical foundation for the use of WH smoothing in constructing experience tables for life insurance risks.

**How to improve the accuracy of smoothing with limited data volume?**

The standard approach applies WH smoothing to crude rate estimates, assuming they are asymptotically normal. However, this assumption often breaks down in practice when data are limited, making the method unreliable in such cases. Following Verrall (1993), we propose
### 1.7.1 How to measure uncertainty in smoothing results?

We propose to measure the uncertainty in smoothing results based on the available observation volume. In a frequentist approach, the estimator associated with smoothing appears to be biased, which makes it challenging to construct relevant confidence intervals for it. However, smoothing can be naturally interpreted in a Bayesian framework, allowing the construction of credible intervals.

### 1.7.2 Which observation and weight vectors to use?

Candidates are proposed for the vectors $y$ and $w$. By adopting the framework of duration models used for experience rating and considering the maximum likelihood estimator of crude rates under the assumption of piecewise constant transition intensity, we present a couple ($y, w$) that satisfies the implicit conditions associated with the use of smoothing in an asymptotic sense.

### 1.7.3 How to improve the accuracy of smoothing with limited data volume?

The previous theoretical framework is based on an asymptotic approximation, which, in many practical situations where the available data volume is limited, can be biased. We propose a generalization of Whittaker-Henderson smoothing based on exact likelihood rather than its normal approximation. This approach requires iterative solving of penalized likelihood equations. By comparing this exact approach with the normal approximation associated with the initial smoothing, we show that this approximation can sometimes come at the cost of substantial bias in the obtained results.

### 1.7.4 How to choose the smoothing parameter(s)?

We then address the choice of the smoothing parameter $\lambda$ to use. This parameter controls the strength of the applied smoothing and has a crucial impact on the obtained results. To select this parameter, we propose maximizing a marginal likelihood using the optimization methods of Brent (1973) and John A. Nelder and Mead (1965), employed in the one-dimensional and two-dimensional case, respectively. These methods are implemented through the `optimize` and `optim` functions in the `stats` library of the statistical programming language R. In the context of the proposed generalization of smoothing, where solving likelihood equations and selecting the smoothing parameter are performed through nested iterations, we compare the precision and computational time of two known strategies called *outer iteration* and *performance iteration*. These strategies differ in the order in which the iterations are nested, and we demonstrate that the latter approach significantly reduces computation time at the cost of negligible loss of precision.
a generalization of WH smoothing that replaces the two-step procedure with the direct maximization of a penalized log-likelihood. Instead of smoothing pre-estimated rates, this method works directly with aggregated event and exposure counts. The estimation is performed iteratively using the PIRLS algorithm. We evaluate both methods on simulated datasets reflecting typical life insurance portfolios. Results show that, in smaller samples, the normal approximation in the traditional method introduces notable bias. This supports the use of the generalized approach—based on penalized log-likelihood—as a more robust alternative when data are limited.

## How to select the smoothing parameters?

We now turn to the crucial choice of the smoothing parameter $\lambda$, which has long been left to actuarial judgment. Giesecke and Center (1981) suggested choosing $\lambda$ so that the variance of the smoothed results matches the average variance of a Chi-square statistic, but uses $n-q$ as degrees of freedom, thus ignoring the reduction in effective model dimension due to penalization. Brooks et al. (1988) minimized the global cross-validation criterion introduced by Wahba (1980), though this can result in severe under-smoothing as noted by Wood (2011).

We instead propose to select $\lambda$ by maximizing a marginal likelihood function, a method first introduced by Patterson and Thompson (1971) and later applied to smoothing parameter selection by Anderssen and Bloomfield (1974). This approach is consistent with the Bayesian framework discussed earlier and performs well in small samples, as shown by Reiss and Todd Ogden (2009). This marginal likelihood function has a closed-form expression and can be maximized numerically. For the proposed generalization of WH smoothing, the marginal likelihood is no longer available in closed form. Instead, we rely on the Laplace Approximation of the Marginal Likelihood (LAML), which can be maximized numerically. As both solving likelihood equations and selecting the optimal smoothing parameter are iterative processes, we explore different ways of nesting these iterations. We compare three nesting strategies combined with three numerical optimization algorithms for maximizing the marginal likelihood or LAML. Simulation results show that all strategies have near-optimal accuracy, with the fastest performance achieved using the outer iteration strategy combined with the Newton algorithm

## How to improve smoothing computational efficiency?

When the number of observations—and thus parameters—is large, the computational cost of WH smoothing becomes a major challenge. This is particularly relevant in actuarial contexts, such as smoothing two-dimensional tables for disability or long-term care modelling. Beyond actuarial applications, WH smoothing is also widely used in economics for long time series, where it is known as the Hodrick-Prescott filter (Hodrick and Prescott 1997). Although fast algorithms have been developed to exploit the structure of the penalization matrix (e.g., Weinert
### 1.7.5 How to improve smoothing performance with a large number of data points?

We propose a practical solution to accelerate smoothing when the number of observations - and thus parameters - is substantial. This situation is common in practice, especially for two-dimensional tables used in the modeling of disability and long-term care risks. To address this issue, an eigendecomposition of the one-dimensional penalization matrices involved in smoothing is used. This decomposition provides a new interpretation of smoothing and opens the door to the use of a rank reduction method to reduce the number of parameters used and the computational time, while incurring minimal loss of precision.

### 1.7.6 How to extrapolate smoothing results?

Finally, we cover the problem of extrapolation in smoothing, which involves solving a modified version of the initial optimization problem where the fidelity criterion remains unchanged and the smoothness criterion is extended to extrapolated values. In the one-dimensional case, a solution may immediately be obtained. However, in the two-dimensional case, as shown by Carballo, Durban, and Lee (2021), constraints must be imposed on the model to prevent the extrapolation step from altering the results obtained during the initial use of the smoothing.

## 2 How to measure uncertainty in smoothing results?

Equation 3 indicates that $E(\hat{y}) = (W + P_{\lambda})^{-1}WE(y) \neq E(y)$ when $\lambda \neq 0$. This implies that penalization introduces a smoothing bias, which prevents the construction of a confidence interval centered around $E(y)$. Therefore, in this section, we turn to a Bayesian approach where smoothing can be interpreted more naturally.

### 2.1 Maximum a posteriori

Let us suppose that $y|\theta \sim N(\theta, W^{-1})$ and $\theta \sim N(0, P_{\lambda}^{-1})$. The Bayes' formula allows us to express the posterior likelihood $f(\theta|y)$ associated with these choices in the following form:

$$ 
\begin{aligned}
f(\theta|y) &= \frac{f(y|\theta)f(\theta)}{f(y)} \\
&\propto f(y|\theta)f(\theta) \\
&\propto \exp\left(-\frac{1}{2}(y-\theta)^T W (y-\theta)\right) \exp\left(-\frac{1}{2}\theta^T P_{\lambda} \theta\right) \\
&\propto \exp\left(-\frac{1}{2}\left[(y-\theta)^T W (y-\theta) + \theta^T P_{\lambda} \theta\right]\right).
\end{aligned}
$$
2007; Cornea-Madeira 2017) they are typically limited to the one-dimensional case and cannot be directly extended to two dimensions.

After briefly outlining the main computational steps of (generalized) WH smoothing-including smoothing parameter selection via marginal likelihood or LAML-and their leading-order costs, we introduce two complementary strategies to reduce the computational burden:

1. Banded matrix exploitation: WH smoothing involves model and penalization matrices with banded structure. Taking advantage of this structure greatly accelerates key computations.

2. Reduced-rank basis via natural parametrization: Building on the work of Demmler and Reinsch (1975), we apply an eigendecomposition to the one-dimensional penalization matrices and drop components associated with the largest eigenvalues, which reduces the problem size. In two dimensions, we further improve efficiency using the Generalized Linear Array Model (GLAM) framework (Currie, Durban, and Eilers 2006) which leverages the rectangular shape of the data.

In the two-dimensional case, we compare these strategies with a cubic P-spline alternative using simulated datasets. Results show that the banded implementation reduces computation time by up to a factor of 25. The reduced-rank approach brings further gains—up to a factor of 250—at the cost of a slight reduction in accuracy. Its performance is comparable to P-spline smoothing with a cubic basis of similar size.

## How to extrapolate smoothing results?

We conclude by addressing how to extrapolate smoothing results. Semi-parametric models like WH and P-splines can extrapolate beyond the observed data—similar to parametric models—but this feature is often overlooked in actuarial practice. The existing literature is limited and mostly focused on mortality forecasting.

Currie, Durban, and Eilers (2004) uses P-splines to fit and forecast mortality rates by treating the extrapolated positions as zero-weight observations (see also Delwarde, Denuit, and Eilers 2007; Currie 2013). While this works well in one dimension, Carballo, Durban, and Lee (2021) showed that it distorts the fit in two dimensions. To fix this, they proposed adding constraints to preserve the values that would result from fitting the observed data alone.

However, their approach to confidence intervals overlooks potential innovation error beyond the observed data, effectively treating the extrapolated process as perfectly smooth. In contrast, we propose an approach that derives credible intervals for extrapolated values, accounting for the underlying variability beyond the observed data range.
The mode of the posterior distribution, $\hat{\theta} = \operatorname{argmax}[f(\theta|\mathbf{y})]$, also known as the maximum a posteriori (MAP) estimate, coincides with the solution $\hat{\mathbf{y}}$ from Equation 2, whose explicit form is given by Equation 3.

## 2.2 Posterior distribution of $\theta|\mathbf{y}$

A second-order Taylor expansion of the log-posterior likelihood around $\hat{\mathbf{y}} = \hat{\theta}$ gives us:

$$ \ln f(\theta|\mathbf{y}) = \ln f(\hat{\theta}|\mathbf{y}) + \left. \frac{\partial \ln f(\theta|\mathbf{y})}{\partial \theta} \right|_{\theta=\hat{\theta}}^T (\theta-\hat{\theta}) + \frac{1}{2}(\theta-\hat{\theta})^T \left. \frac{\partial^2 \ln f(\theta|\mathbf{y})}{\partial \theta \partial \theta^T} \right|_{\theta=\hat{\theta}} (\theta-\hat{\theta}) \quad (4) $$

$$ \text{where} \quad \left. \frac{\partial \ln f(\theta|\mathbf{y})}{\partial \theta} \right|_{\theta=\hat{\theta}} = 0 \quad \text{and} \quad \left. \frac{\partial^2 \ln f(\theta|\mathbf{y})}{\partial \theta \partial \theta^T} \right|_{\theta=\hat{\theta}} = -(W + P_{\lambda}). $$

Note that this last derivative no longer depends on $\theta$, and higher-order derivatives beyond 2 are all zero. The Taylor expansion allows for an exact computation of $\ln f(\theta|\mathbf{y})$. By substituting the result back into Equation 4, we obtain:

$$ \begin{aligned} f(\theta|\mathbf{y}) &\propto \exp \left[ \ln f(\hat{\theta}|\mathbf{y}) - \frac{1}{2}(\theta-\hat{\theta})^T (W + P_{\lambda}) (\theta-\hat{\theta}) \right] \\ &\propto \exp \left[ -\frac{1}{2}(\theta-\hat{\theta})^T (W + P_{\lambda}) (\theta-\hat{\theta}) \right] \end{aligned} $$

which can immediately be recognized as the density of the $N(\hat{\theta}, (W + P_{\lambda})^{-1})$ distribution.

## 2.3 Consequence for Whittaker-Henderson smoothing

The assumption $\theta \sim N(0, P_{\lambda}^{-1})$ corresponds to a simple Bayesian formalization of the smoothness criterion. It reflects a (improper) prior belief of the modeler about the underlying distribution of the observation vector $\mathbf{y}$.

The use of Whittaker-Henderson smoothing in the Bayesian framework and the construction of credible intervals are conditioned on the validity of the assumption $\mathbf{y}|\theta \sim N(\theta, W^{-1})$. The components of the observation vector should be independent and have known variances. The weight vector $\mathbf{w}$ used should contain the inverses of these variances. Under these assumptions, credible intervals at $100(1 - \alpha)\%$ for smoothing can be constructed and take the form:

$$ \mathbb{E}(\mathbf{y})|\mathbf{y} \in \left[ \hat{\mathbf{y}} \pm \Phi^{-1}\left(1 - \frac{\alpha}{2}\right) \sqrt{\operatorname{diag}\{(W + P_{\lambda})^{-1}\}} \right] \quad (5) $$
## 2 How to measure uncertainty in smoothing results?

The explicit solution given by Equation 3 indicates that $E(\hat{y}) = (W + P_\lambda)^{-1}WE(y) \neq E(y)$ when $\lambda \neq 0$. This implies that penalization introduces a smoothing bias, which prevents the construction of confidence intervals for finite samples centred on $E(y)$. Therefore, in this section, we turn to a Bayesian framework where smoothing can be interpreted more naturally.

### 2.1 Maximum a posteriori estimate

Suppose that $y | \theta \sim N(\theta, \sigma^2 W^{-1})$ and $\theta \sim N(0, \sigma^2 P_\lambda^{-1})$ for some $\sigma > 0$. The Bayes formula allows us to express the posterior likelihood $f(\theta | y)$ associated with these choices in the following form:

$$f(\theta | y) \propto f(y | \theta)f(\theta) \propto \exp \left( -\frac{1}{2\sigma^2} \left[ (y - \theta)^T W (y - \theta) + \theta^T P_\lambda \theta \right] \right).$$

Hence the mode of the posterior distribution, $\hat{\theta} = \operatorname{argmax}[f(\theta | y)]$, also known as the maximum a posteriori (MAP) estimate, coincides with the solution $\hat{y}$ from Equation 2, whose explicit form is given by Equation 3.

### 2.2 Posterior distribution of $\theta$ | y

A second-order Taylor expansion of the log-posterior likelihood around $\hat{y} = \hat{\theta}$ gives us:

$$\ln f(\theta | y) = \ln f(\hat{\theta} | y) + \left. \frac{\partial \ln f(\theta | y)}{\partial \theta} \right|_{\theta=\hat{\theta}}^T (\theta - \hat{\theta}) + \frac{1}{2} (\theta - \hat{\theta})^T \frac{\partial^2 \ln f(\theta | y)}{\partial \theta \partial \theta^T} \bigg|_{\theta=\hat{\theta}} (\theta - \hat{\theta}) \quad (4)$$

$$\text{where } \left. \frac{\partial \ln f(\theta | y)}{\partial \theta} \right|_{\theta=\hat{\theta}} = 0 \quad \text{and} \quad \left. \frac{\partial^2 \ln f(\theta | y)}{\partial \theta \partial \theta^T} \right|_{\theta=\hat{\theta}} = -\frac{1}{\sigma^2}(W + P_{\lambda}).$$

As this last derivative no longer depends on $\theta$, higher-order derivatives are all zero. The Taylor expansion allows for an exact computation of $\ln f(\theta | y)$. Substituting the result back into Equation 4 yields:

$$f(\theta | y) \propto \exp \left[ \ln f(\hat{\theta} | y) - \frac{1}{2\sigma^2} (\theta - \hat{\theta})^T (W + P_{\lambda}) (\theta - \hat{\theta}) \right] \\ \propto \exp \left[ -\frac{1}{2\sigma^2} (\theta - \hat{\theta})^T (W + P_{\lambda}) (\theta - \hat{\theta}) \right]$$

which can immediately be recognized as the density of the $N(\hat{\theta}, \sigma^2(W + P_\lambda)^{-1})$ distribution.
with probability $1 - \frac{\alpha}{2}$ where $\hat{y} = (W + P_\lambda)^{-1}W y$ and $\Phi$ denotes the cumulative distribution function of the standard normal distribution. According to Marra and Wood (2012), these credible intervals provide satisfactory coverage of the corresponding frequentist confidence intervals and can therefore be used as practical substitutes.

# 3 Which observation and weight vectors to use?

Section 2 highlighted the need to apply the Whittaker-Henderson smoothing to independent observation vectors $y$ and weight vectors $w$ corresponding to the inverses of the variances of the components of $y$ in order to obtain a measure of uncertainty in the results. In this section, we propose, within the framework of duration models used for constructing experience tables, vectors $y$ and $w$ that satisfy these conditions.

## 3.1 Duration models framework: one-dimensional case

Consider the observation of $m$ individuals in a longitudinal study subject to left truncation and right censoring phenomena. Suppose we want to estimate a single distribution that depends on only one continuous explanatory variable, denoted by $x$. For illustration purposes, we can consider it as a mortality distribution with the explanatory variable of interest $x$ representing age. Such a distribution is fully characterized by one of the following quantities:

* the cumulative distribution function $F(x)$ or its complement, the survival function $S(x) = 1 - F(x)$,

* the associated probability density function $f(x) = -\frac{d}{dx}S(x)$,

* the instantaneous hazard function $\mu(x) = -\frac{d}{dx}\ln S(x)$.

Suppose that the considered distribution depends on a vector of parameters $\theta$ that we want to estimate using maximum likelihood. The likelihood associated with the observation of the individuals can be written as follows:

$$ \mathcal{L}(\theta) = \prod_{i=1}^{m} \left[ \frac{f(x_i + t_i, \theta)}{S(x_i, \theta)} \right]^{\delta_i} \left[ \frac{S(x_i + t_i, \theta)}{S(x_i, \theta)} \right]^{1-\delta_i} \quad (6) $$

Where $x_i$ represents the age at the start of observation, $t_i$ represents the observation duration, i.e., the time elapsed between the start and end dates of observation, and $\delta_i$ is the indicator of event observation, which takes the value 1 if the event of interest is observed and 0 if the observation is censored. We will not go into the details of calculating these three quantities, which should take into account the individual-specific information such as the subscription date, redemption date if applicable, as well as the global characteristics of the product, such as the presence of a deductible or medical selection phenomenon, and the choice of a restricted observation period due to data quality issues or delays in the reporting of observations. These
## 2.3 Consequence for the WH smoothing

The prior $\theta \sim \mathcal{N}(0, \sigma^2 P_{\lambda}^{-1})$ provides a Bayesian interpretation of the smoothness penalty, expressing an (improper) prior belief about the structure of $\mathbf{y}$.

This Bayesian framework and the resulting credible intervals rely on the assumption that $\mathbf{y} | \theta \sim \mathcal{N}(\theta, \sigma^2 W^{-1})$, meaning that the components of $\mathbf{y}$ are independent with known variances (up to a constant $\sigma^2$). The weight vector $\mathbf{w}$ must then be proportional to the inverse variances, not chosen empirically. If $\sigma^2$ is known, $100(1 - \alpha)\%$ credible intervals take the form:

$$ \mathbb{E}(\mathbf{y}) | \mathbf{y} \in \left[ \hat{\mathbf{y}} \pm \Phi^{-1} \left(1 - \frac{\alpha}{2}\right) \sqrt{\sigma^2 \mathrm{diag}\{(W + P_{\lambda})^{-1}\}} \right] \quad (5) $$

where $\hat{\mathbf{y}} = (W + P_{\lambda})^{-1} W \mathbf{y}$ and $\Phi$ is the cumulative distribution function for the standard normal distribution. According to Marra and Wood (2012), such intervals have good Frequentist coverage.

If $\sigma^2$ is unknown, it can be estimated as:

$$ \hat{\sigma}^2 = \frac{(\mathbf{y} - \hat{\mathbf{y}})^T W (\mathbf{y} - \hat{\mathbf{y}})}{n - \mathrm{tr}(H)} \quad \text{where} \quad H = (W + P_{\lambda})^{-1} W. $$

In that case, $\sigma^2$ is replaced by $\hat{\sigma}^2$ and the normal distribution in Equation 5 by the Student t - distribution with $n - \mathrm{tr}(H)$ degrees of freedom.

# 3 Which observation and weight vectors to use?

Section 2 highlighted that Whittaker-Henderson smoothing may be interpreted in a robust statistical framework when applied to a vector $\mathbf{y}$ of independent, normally distributed observations with known variances, and a weight vector $\mathbf{w}$ proportional to the inverses of those variances. In this section, we propose, within the framework of duration models used for constructing experience tables for life insurance risks, vectors $\mathbf{y}$ and $\mathbf{w}$ that satisfy these conditions.

## 3.1 Survival analysis framework

We consider a longitudinal follow-up of $m$ individuals, subject to left truncation and non-informative right censoring, and aim to estimate a distribution governed by a continuous explanatory variable $x$ (e.g., age). Let $\mu$ denote the hazard function, also known as the force of mortality in the study of the death risk. Under standard survival analysis assumptions, the log-likelihood takes the following continuous-time form:

$$ \ell(\theta) = \sum_{i=1}^{m} \left[ \delta_i \ln \mu(x_i + t_i, \theta) - \int_{u=0}^{t_i} \mu(x_i + u, \theta) du \right]. \quad (6) $$
factors typically lead to a narrower observation period than the actual presence period of individuals in the portfolio.

The various quantities introduced above are related by the following relationships:

$$S(x) = \exp \left( \int_{u=0}^{x} \mu(u) du \right) \quad \text{and} \quad f(x) = \mu(x)S(x).$$

The maximization of the likelihood in Equation 6 is equivalent to the maximization of the associated log-likelihood, which can be rewritten using only the instantaneous hazard function (also known as the mortality rate in the case of the death risk):

$$\ell(\boldsymbol{\theta}) = \sum_{i=1}^{m} \left[ \delta_i \ln \mu(x_i + t_i, \boldsymbol{\theta}) - \int_{u=0}^{t_i} \mu(x_i + u, \boldsymbol{\theta}) du \right] \quad (7)$$

To discretize the problem, let us assume that the mortality rate is piecewise constant over one-year intervals between two integer ages. Formally, we have $\mu(x + \epsilon) = \mu(x)$ for all $x \in \mathbb{N}$ and $\epsilon \in [0, 1[$. Furthermore, if $\mathbf{1}$ denotes the indicator function, then for any $0 \le a < x_{\max}$, we have $\sum_{x=x_{\min}}^{x_{\max}} \mathbf{1}(x \le a < x + 1) = 1$, where $x_{\min} = \min(\mathbf{x})$ and $x_{\max} = \max(\mathbf{x})$. Equation 7 can be rewritten as:

$$\ell(\boldsymbol{\theta}) = \sum_{i=1}^{m} \left[ \sum_{x=x_{\min}}^{x_{\max}} \delta_i \mathbf{1}(x \le x_i + t_i < x + 1) \ln \mu(x_i + t_i, \boldsymbol{\theta}) - \int_{u=0}^{t_i} \sum_{x=x_{\min}}^{x_{\max}} \mathbf{1}(x \le x_i + u < x + 1) \mu(x_i + u, \boldsymbol{\theta}) du \right].$$

The assumption of piecewise constant mortality rate implies that:

$$\mathbf{1}(x \le x_i + t_i < x + 1) \ln \mu(x_i + t_i, \boldsymbol{\theta}) = \mathbf{1}(x \le x_i + t_i < x + 1) \ln \mu(x, \boldsymbol{\theta}) \quad \text{and} \\ \mathbf{1}(x \le x_i + u < x + 1) \mu(x_i + u, \boldsymbol{\theta}) = \mathbf{1}(x \le x_i + u < x + 1) \ln \mu(x, \boldsymbol{\theta}).$$

It is then possible to interchange the two summations to obtain the following expressions:
Here $x_i$ is the age at the start of observation, $t_i$ is the follow-up duration for individual $i$ and $\delta_i$ is an event indicator: 1 if the event is observed and 0 if censored.

Although model estimation can be based on direct maximization of Equation 6, this approach scales poorly with large $m$ and generally requires numerical integration-except in simple parametric cases. We instead adopt a discrete approximation by assuming the hazard rate is piecewise constant over one-year intervals:

$$ \mu(x + \epsilon) = \mu(x) \quad \text{for all } x \in \mathbb{N}, \epsilon \in [0, 1[. $$

Under this assumption, the log-likelihood simplifies to a sum over discrete ages:

$$ \ell(\boldsymbol{\theta}) = \sum_{x=x_{\min}}^{x_{\max}} \ln \mu(x, \boldsymbol{\theta})d(x) - \mu(x, \boldsymbol{\theta})e_c(x). \quad (7) $$

Here $d(x)$ is the number of observed events at age $x$ and $e_c(x)$ is the central exposure to risk, i.e., the total duration individuals are observed at age $x$.

This discretization, first introduced by Hoem (1971), is widely used in actuarial science. Its advantages are underlined for example in Gschlössl, Schoenmaekers, and Denuit (2011). It extends naturally to the two-dimensional case by assuming $\mu(x + \epsilon, z + \xi) = \mu(x, z)$ and summing over $(x, z)$ pairs.

Details on the derivation of Equations 6 and 7, along with the computation of central exposures and event counts, are provided in Section A of the appendices.

## 3.2 Likelihood equations

Assuming one parameter per observation and using the exponential link $\mu(\boldsymbol{\theta}) = \exp(\boldsymbol{\theta})$, we recover the crude rates estimator, which models each age (or age pair) independently. The exponential link ensures positive hazard rates. The log-likelihood, in both one- and two-dimensional cases, takes the vectorized form:

$$ \ell(\boldsymbol{\theta}) = \boldsymbol{\theta}^T \mathbf{d} - \exp(\boldsymbol{\theta})^T \mathbf{e}_c \quad (8) $$

where $\mathbf{d}$ and $\mathbf{e}_c$ are the vectors of observed deaths and central exposures.

The derivatives of this likelihood are:

$$ \frac{\partial \ell}{\partial \boldsymbol{\theta}} = \mathbf{d} - \exp(\boldsymbol{\theta}) \odot \mathbf{e}_c \quad \text{and} \quad \frac{\partial^2 \ell}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^T} = -\text{Diag}(\exp(\boldsymbol{\theta}) \odot \mathbf{e}_c). \quad (9) $$

These equations correspond to those of a Poisson GLM (Nelder and Wedderburn 1972) with mean $\mu(\boldsymbol{\theta}) \odot \mathbf{e}_c$, although derived under different assumptions.
$$ \ell(\boldsymbol{\theta}) = \sum_{x=x_{\min}}^{x_{\max}} [\ln \mu(x, \boldsymbol{\theta})d(x) - \mu(x, \boldsymbol{\theta})e_c(x)] \quad \text{where} $$

$$ d(x) = \sum_{i=1}^{m} \delta_i \mathbf{1}(x \le x_i + t_i < x + 1) \quad \text{and} $$

$$ e_c(x) = \sum_{i=1}^{m} \int_{u=0}^{t_i} \mathbf{1}(x \le x_i + u < x + 1) du = \sum_{i=1}^{m} [\min(t_i, x - x_i + 1) - \max(0, x - x_i)]^+ $$

by denoting $a^+ = \max(a, 0)$, where $d(x)$ and $e_c(x)$ correspond to the number of observed deaths between ages $x$ and $x+1$ and the sum of observation durations of individuals between these ages, respectively (the latter quantity is also known as central exposure to risk).

## 3.2 Extension to the two-dimensional case

The extension of the proposed approach to the two-dimensional framework requires only minor adjustments to the previous reasoning. Let $z_{\min} = \min(\mathbf{z})$ and $z_{\max} = \max(\mathbf{z})$. The piecewise constant assumption for the mortality rate needs to be extended to the second dimension. Formally, we now assume that $\mu(x+\epsilon, z+\xi) = \mu(x, z)$ for all pairs $x, z \in \mathbb{N}$ and $\epsilon, \xi \in [0, 1[$. The sums involving the variable $x$ are then replaced by double sums considering all combinations of $x$ and $z$. The log-likelihood is given by:

$$ \ell(\boldsymbol{\theta}) = \sum_{x=x_{\min}}^{x_{\max}} \sum_{z=z_{\min}}^{z_{\max}} [\ln \mu(x, z, \boldsymbol{\theta})d(x, z) - \mu(x, z, \boldsymbol{\theta})e_c(x, z)] \quad \text{where} $$

$$ d(x, z) = \sum_{i=1}^{m} \delta_i \mathbf{1}(x \le x_i + t_i < x + 1) \mathbf{1}(z \le z_i + t_i < z + 1) \quad \text{and} $$

$$ \mathbf{e}_c(x, z) = \sum_{i=1}^{m} \int_{u=0}^{t_i} \mathbf{1}(x \le x_i + u < x + 1) \mathbf{1}(z \le z_i + u < z + 1) du $$

$$ = \sum_{i=1}^{m} [\min(t_i, x + 1 - x_i, z + 1 - z_i) - \max(0, x - x_i, z - z_i)]^+ $$

## 3.3 Likelihood equations

The choice $\mu(\boldsymbol{\theta}) = \exp(\boldsymbol{\theta})$, which corresponds to considering one parameter per observation, allows us to relate to the Whittaker-Henderson smoothing. Using the exponential function ensures positive values for the estimated mortality rate. The expressions of the likelihood in the one-dimensional or two-dimensional case can then be written in a common vectorized form:
The model admits the closed-form solution $\hat{\theta} = \ln(\mathbf{d}/\mathbf{e}_c)$. Under standard regularity conditions, the maximum likelihood estimator satisfies $\hat{\theta} \sim \mathcal{N}(\theta, W_{\hat{\theta}}^{-1})$, with $W_{\hat{\theta}} = \text{Diag}(\mathbf{d})$.

Notably, this asymptotic approximation depends on the number of individuals $m$ and not the dimension $n$ of the aggregated vectors.

### 3.3 Consequence for the WH smoothing

We conclude that, under the duration model framework and using crude rates, the log-estimate $\ln(\mathbf{d}/\mathbf{e}_c)$ is asymptotically normal:

$$ \ln(\mathbf{d}/\mathbf{e}_c) \sim \mathcal{N}(\ln \boldsymbol{\mu}, W^{-1}) \quad \text{with} \quad W = \text{Diag}(\mathbf{d}). $$

This justifies applying Whittaker-Henderson smoothing to the observation vector $\mathbf{y} = \ln(\mathbf{d}/\mathbf{e}_c)$ with weight vector $\mathbf{w} = \mathbf{d}$. Using results from Section 2, and $\sigma^2 = 1$, the credible intervals for $\ln \boldsymbol{\mu}$ are:

$$ \ln \boldsymbol{\mu} | \mathbf{d}, \mathbf{e}_c \in \left[ \hat{\boldsymbol{\theta}} \pm \Phi^{-1} (1 - \alpha/2) \sqrt{\text{diag} \{ (\text{Diag}(\mathbf{d}) + P_{\lambda})^{-1} \}} \right] $$

with $\hat{\boldsymbol{\theta}} = (W + P_{\lambda})^{-1}W(\ln \mathbf{d} - \ln \mathbf{e}_c)$. Credible intervals for $\boldsymbol{\mu}$ itself are then obtained by exponentiating the bounds.

## 4 How to improve the accuracy of smoothing with limited data volume?

### 4.1 Generalized Whittaker-Henderson smoothing

The approach described in Section 3.2 assumes that the crude rates estimator is asymptotically normal, justifying the application of WH smoothing to its logarithm. However, with limited data, this approximation may introduce significant bias. We therefore propose an alternative based directly on the exact likelihood in Equation 8. Applying the Bayesian framework from Section 2 and assuming $\theta \sim \mathcal{N}(0, P_{\lambda}^{-1})$, Bayes' theorem gives:

$$ f(\boldsymbol{\theta} | \mathbf{d}, \mathbf{e}_c) \propto f(\mathbf{d}, \mathbf{e}_c | \boldsymbol{\theta}) f(\boldsymbol{\theta}) \propto \exp \left[ \ell(\boldsymbol{\theta}) - \frac{1}{2} \boldsymbol{\theta}^T P_{\lambda} \boldsymbol{\theta} \right]. $$

We define the penalized log-likelihood as $\ell_P(\boldsymbol{\theta}) = \ell(\boldsymbol{\theta}) - \boldsymbol{\theta}^T P_{\lambda} \boldsymbol{\theta}/2$. The maximum a posteriori estimate is the maximizer of $\ell_P$.

Using a second-order Taylor expansion of the posterior log-likelihood around $\hat{\boldsymbol{\theta}}$ leads to the Laplace approximation:
$$ \ell(\theta) = \theta^T d - \exp(\theta)^T e_c \qquad (8) $$

where d and e_c represent the vectors of observed deaths and exposures to risk, respectively. The derivatives of the likelihood function for this model are given by:

$$ \frac{\partial \ell}{\partial \theta} = [d - \exp(\theta) \odot e_c] \quad \text{and} \quad \frac{\partial^2 \ell}{\partial \theta \partial \theta^T} = -\text{Diag}(\exp(\theta) \odot e_c). \qquad (9) $$

Note that these likelihood equations are exactly what we would obtain by assuming that the observed numbers of deaths, conditional on the observed exposures to risk e_c, follow Poisson distributions with parameters μ(θ) ⊙ e_c. The model presented here has many similarities with a Poisson GLM (John Ashworth Nelder and Wedderburn 1972). However, the initial assumptions are not the same for both models.

These likelihood equations have an explicit solution given by $\hat{\theta} = \ln(d/e_c)$. This model, which treats each age independently, is known as the crude rates estimator. The properties of the maximum likelihood estimator imply that asymptotically $\hat{\theta} \sim N(\theta, W^{-1})$, where $W_{\hat{\theta}}$ is a diagonal matrix with elements $\exp(\hat{\theta}) \odot e_c = (d/e_c) \odot e_c = d$. It should be noted that the asymptotic nature and the validity of this approximation are related to the number of individuals m in the portfolio and not the size n of the vectors d and e_c.

Thus, we have shown that in the framework of duration models, using the crude rates estimator, asymptotically $\ln(d/e_c) \sim N(\ln \mu, W^{-1})$, where $W = \text{Diag}(d)$. This justifies applying the Whittaker-Henderson smoothing to the observation vector $y = \ln(d/e_c)$ and weight vector $w = d$. According to the results from Section 2, we obtain the following associated credible intervals:

$$ \ln \mu | d, e_c \in [\hat{\theta} \pm \Phi^{-1} \left(1 - \frac{\alpha}{2}\right) \sqrt{\text{diag} \{ (\text{Diag}(d) + P_{\lambda})^{-1} \}}] $$

where $\hat{\theta} = (\text{Diag}(d) + P_{\lambda})^{-1}\text{Diag}(d)[\ln(d) - \ln e_c]$. Results on μ can be obtained directly by exponentiating the above expression.

# 4 How to improve the accuracy of smoothing with limited data volume?

## 4.1 Generalized Whittaker-Henderson smoothing

The approach described in Section 3.3 relies on the asymptotic properties of the maximum likelihood estimator, which provide a theoretical framework for applying smoothing using the
$$f(\theta | \mathbf{d}, \mathbf{e}_c) \approx \mathcal{N}(\hat{\theta}, (W_{\hat{\theta}} + P_{\lambda})^{-1}) \quad (10)$$

where $W_{\hat{\theta}} = \text{Diag}(\exp(\hat{\theta}) \odot \mathbf{e}_c)$. Unlike the normal case studied in Section 2, the higher-order derivatives of the posterior log-likelihood are not zero, and Equation 10 only provides an approximation of the posterior log-likelihood, which yields asymptotic credible intervals:

$$\ln \mu | \mathbf{d}, \mathbf{e}_c \in [\hat{\theta} \pm \Phi^{-1}(1 - \alpha/2) \sqrt{\text{diag} \{(W_{\hat{\theta}} + P_{\lambda})^{-1}\}}].$$

Unlike the closed-form estimator in Equation 9, no analytical solution for $\hat{\theta}$ exists here. We solve numerically using Newton's algorithm, which iteratively updates:

$$\theta_{k+1} = \theta_k + (W_k + P_{\lambda})^{-1}(\mathbf{d} - \exp(\theta_k) \odot \mathbf{e}_c - P_{\lambda}\theta_k)$$

with $W_k = \text{Diag}(\exp(\theta_k) \odot \mathbf{e}_c)$. The update can be rewritten as:

$$\theta_{k+1} = (W_k + P_{\lambda})^{-1} W_k \mathbf{z}_k \quad \text{where} \quad \mathbf{z}_k = \theta_k + W_k^{-1} [\mathbf{d} - \exp(\theta_k) \odot \mathbf{e}_c].$$

Initializing with the crude rates estimator $\theta_0 = \ln(\mathbf{d}/\mathbf{e}_c)$ implies $W_0 = \text{Diag}(\mathbf{d})$ and $z_0 = \ln(\mathbf{d}/\mathbf{e}_c)$, so the first iteration recovers the classical WH smoothing result.

Subsequent iterations refine the observation and weight vectors. This process can thus be interpreted as an iterative generalization of WH smoothing, akin to how generalized linear models extend linear models.

We refer to this method as generalized Whittaker-Henderson smoothing. The iterative estimation algorithm described above corresponds to the Penalized Iteratively Reweighted Least Squares (PIRLS) algorithm, widely used for fitting generalized additive models.

This framework naturally extends to other exponential family distributions, such as the binomial case suggested in Verrall (1993), by adapting the likelihood, link function, weight matrix, and working vector. However, we advocate for the Poisson-like likelihood of Equation 8, which offers several advantages: it generalizes to competing risks, supports multiplicative covariate effects via the log link, and allows the use of an external reference table as a multiplicative offset.

## 4.2 Impact of the normal approximation in the original smoothing

As discussed in Section 3, classical Whittaker-Henderson smoothing can be viewed as an approximation to a penalized likelihood maximization, relying on a crude rate estimator assumed to be asymptotically normal. To assess the practical consequences of this approximation, we conduct an empirical comparison based on six simulated datasets reflecting the typical structure and volume of real insurance portfolios:
crude rates estimator. However, the validity of these asymptotic properties in practice, where we have a limited number of observations, may be questioned.

Another approach is to apply the Bayesian reasoning presented in Section 2 directly to the likelihood of Equation 8. Let us assume again that $\theta \sim \mathcal{N}(0, P_{\lambda}^{-1})$ and write, using Bayes' theorem:

$$f(\theta|\mathbf{d}, \mathbf{e}_c) = \frac{f(\mathbf{d}, \mathbf{e}_c|\theta)f(\theta)}{f(\mathbf{d}, \mathbf{e}_c)} \propto f(\mathbf{d}, \mathbf{e}_c|\theta)f(\theta) \propto \exp\left[\ell(\theta) - \frac{1}{2}\theta^T P_{\lambda} \theta\right].$$

The quantity $\ell_P(\theta) = \ell(\theta) - \frac{1}{2}(\theta)^T P_{\lambda} (\theta)$ will be referred to as the penalized likelihood. The maximum a posteriori $\hat{\theta} = \underset{\theta}{\operatorname{argmax}} f(\theta|\mathbf{d}, \mathbf{e}_c)$ corresponds, once again, to the maximum of the penalized likelihood $\tilde{\theta} = \underset{\theta}{\operatorname{argmax}} \ell_P(\theta)$.

A second-order Taylor expansion of the posterior log-likelihood around $\hat{\theta}$ yields:

$$\ln f(\theta|\mathbf{d}, \mathbf{e}_c) \approx \ln f(\hat{\theta}|\mathbf{d}, \mathbf{e}_c) + \left. \frac{\partial \ln f(\theta|\mathbf{d}, \mathbf{e}_c)}{\partial \theta} \right|_{\theta=\hat{\theta}}^T (\theta-\hat{\theta}) + \frac{1}{2} (\theta-\hat{\theta})^T \left. \frac{\partial^2 \ln f(\theta|\mathbf{d}, \mathbf{e}_c)}{\partial \theta \partial \theta^T} \right|_{\theta=\hat{\theta}} (\theta-\hat{\theta}) \quad (10)$$

$$\text{avec ici } \left. \frac{\partial \ln f(\theta|\mathbf{d}, \mathbf{e}_c)}{\partial \theta} \right|_{\theta=\hat{\theta}} = 0 \quad \text{and} \quad \left. \frac{\partial^2 \ln f(\theta|\mathbf{d}, \mathbf{e}_c)}{\partial \theta \partial \theta^T} \right|_{\theta=\hat{\theta}} = -(W_{\hat{\theta}} + P_{\lambda}).$$

where $W_{\hat{\theta}} = \operatorname{Diag}(\exp(\hat{\theta}) \odot \mathbf{e}_c)$.

Unlike the normal case studied in Section 2, the higher-order derivatives of the posterior log-likelihood are not zero, and the Equation 10 represents an approximation of the posterior log-likelihood known as the Laplace approximation, which becomes more accurate as the number of observations increases. Asymptotically:

$$f(\theta|\mathbf{d}, \mathbf{e}_c) \propto \exp\left[\ell_P(\hat{\theta}) - \frac{1}{2}(\theta - \hat{\theta})^T(W_{\hat{\theta}} + P_{\lambda})(\theta - \hat{\theta})\right] \\ \propto \exp\left[-\frac{1}{2}(\theta - \hat{\theta})^T(W_{\hat{\theta}} + P_{\lambda})(\theta - \hat{\theta})\right]$$

and the posterior distribution is thus equivalent to $\mathcal{N}(\hat{\theta}, (W_{\hat{\theta}} + P_{\lambda})^{-1})$. Since $\theta = \ln \mu$, this result allows for the construction of asymptotic credible intervals at $100(1 - \alpha)$ for $\ln \mu|\mathbf{d}, \mathbf{e}_c$, in the form:

$$\ln \mu|\mathbf{d}, \mathbf{e}_c \in \left[ \hat{\theta} \pm \Phi^{-1}\left(1 - \frac{\alpha}{2}\right) \sqrt{\operatorname{diag}\{(W_{\hat{\theta}} + P_{\lambda})^{-1}\}} \right].$$
<table><caption>Table 1: Key figures associated with the 6 simulated datasets</caption><thead><tr><th>Portfolio type</th><th>Dimensions</th><th>Head count</th><th>Exposure count</th><th>Death count</th></tr></thead><tbody><tr><td>annuity</td><td>45</td><td>20,000</td><td>136,524</td><td>1,722</td></tr><tr><td>annuity</td><td>45</td><td>100,000</td><td>679,728</td><td>8,452</td></tr><tr><td>annuity</td><td>45</td><td>500,000</td><td>3,405,892</td><td>42,499</td></tr><tr><td>LTC</td><td>30 x 15</td><td>20,000</td><td>8,115</td><td>1,888</td></tr><tr><td>LTC</td><td>30 x 15</td><td>100,000</td><td>40,004</td><td>9,281</td></tr><tr><td>LTC</td><td>30 x 15</td><td>500,000</td><td>202,666</td><td>47,358</td></tr></tbody></table>

* The first three datasets simulate annuity portfolios with 20,000, 100,000, and 500,000 policyholders. The sole covariate is age, ranging from 50 to 95.

* The next three mimic long-term care (LTC) portfolios of the same sizes. Modelling of LTC typically relies on the illness-death model (Fix and Neyman 1951; Clifford 1977). To get a two-dimensional illustration we focus on the transition between the disabled and dead states (the two other transitions would provide additional one-dimensional examples). Two covariates are used: age (70–100) and duration in LTC (0–15 years).

Each dataset consists of individual-level longitudinal data, from which we derive event counts $d$ and exposures $e_c$, aggregated by age $x$ (for annuities) of by $(x, z)$ pairs (for LTC). All datasets within each group share the same underlying structure and differ only in size. Key dataset statistics are provided in Table 1 and additional details about how those datasets were generated are provided in Section B of the appendices.

We apply two methods:

1. Original WH smoothing using $y = \ln(d/e_c)$ and weights $w = d$ as in Section 3.

2. Generalized WH smoothing, using the likelihood formulation of Section 4.

Both methods use the same smoothing parameter(s) $\lambda$, to ensure that prior assumptions on $\theta = \ln \mu$ are held constant. We fix the penalty order at $q = 2$, corresponding to second-order differences.

As both estimators target $\theta$, we compare them using the following relative error metric:

$$ \Delta(\theta) = \frac{\ell_P(\hat{\theta}_{\text{ML}}) - \ell_P(\theta)}{\ell_P(\hat{\theta}_{\text{ML}}) - \ell_P(\hat{\theta}_{\infty})} \quad (11) $$

Here $\hat{\theta}_{\infty}$ maximizes the penalized likelihood, while $\hat{\theta}_{\infty}$ corresponds to the solution with $\lambda \to \infty$, which we later show to be the degree-(q - 1) polynomial that maximizes the likelihood. By construction:

$$ \Delta(\hat{\theta}_{\text{ML}}) = 0, \quad \Delta(\hat{\theta}_{\infty}) = 1, \quad \text{and} \quad \Delta(\theta) \ge 0. $$
Unlike Equation 9, it is not possible to determine the maximum of the penalized likelihood $\hat{\theta}$ explicitly. However, the Newton algorithm allows for a numerical solution of the likelihood equations by constructing a sequence of estimators $(\theta_k)_{k \ge 0}$ that converges to $\hat{\theta} = \underset{\theta}{\operatorname{argmax}} \ell_P(\theta)$.

These estimators are recursively defined as:

$$ \begin{aligned} \theta_{k+1} &= \theta_k - \left( \left. \frac{\partial^2 \ell_P}{\partial \theta^2} \right|_{\theta=\theta_k} \right)^{-1} \left. \frac{\partial \ell_P}{\partial \theta} \right|_{\theta=\theta_k} \\ &= \theta_k + (W_k + P_\lambda)^{-1} (\mathbf{d} - \exp(\theta_k) \odot \mathbf{e}_c - P_\lambda \theta_k)] \\ &= (W_k + P_\lambda)^{-1} W_k \mathbf{z}_k \end{aligned} $$

by denoting $W_k = \operatorname{Diag}(\exp(\theta_k) \odot \mathbf{e}_c)$ and $\mathbf{z}_k = \theta_k + W_k^{-1} [\mathbf{d} - \exp(\theta_k) \odot \mathbf{e}_c]$. An interesting initialization choice for the algorithm starts with the crude rates estimator, setting $\theta_0 = \ln(\mathbf{d}/\mathbf{e}_c)$, from which we derive $W_0 = \operatorname{Diag}(\mathbf{d})$ and $z_0 = \ln(\mathbf{d}/\mathbf{e}_c)$.

With this choice $\theta_1 = (W_0 + P_\lambda)^{-1} W_0 z_0$ corresponds to the solution of the original Whittaker-Henderson smoothing. Subsequent iterations can also be seen as successive applications of the original smoothing to pseudo-observation vectors $\mathbf{z}_k$ and associated weight vectors $w_k = \operatorname{diag}(W_k)$, adjusted at each iteration. The proposed approach can thus be considered as an iterative variant of the original Whittaker-Henderson smoothing, where the choice of observation and weight vectors is refined at each step. The connection between the approach presented in this section and the original smoothing is of the same nature as that between a linear model and its generalized version. Therefore, we will refer to this approach as generalized Whittaker-Henderson smoothing throughout the rest of the paper. Algorithm 1 provides an implementation of this new version of the smoothing.

## 4.2 Impact of the normal approximation of the original Whittaker-Henderson smoothing

We have shown that the original Whittaker-Henderson smoothing can be interpreted as a normal approximation of a penalized likelihood maximization problem. To assess the impact of this approximation in practice, we rely on six simulated datasets that capture the characteristics of real biometric risks. These datasets include:

* In the one-dimensional case, three annuity portfolios with 20,000, 100,000, and 500,000 policyholders, respectively.

* In the two-dimensional case, three annuity portfolios for policyholders in a long-term care situation, with 1,000, 5,000, and 25,000 dependent policyholders, respectively.
<table><caption>Table 2: Impact of the approximation from the original WH smoothing on the 6 simulated datasets</caption><thead><tr><td>Portfolio type</td><td>Head count</td><td>Relative Error</td><td>SMR</td></tr></thead><tbody><tr><td>annuity</td><td>20,000</td><td>1,91%</td><td>99,19%</td></tr><tr><td>annuity</td><td>100,000</td><td>0,02%</td><td>99,89%</td></tr><tr><td>annuity</td><td>500,000</td><td>0,00%</td><td>99,99%</td></tr><tr><td>LTC</td><td>20,000</td><td>93,27%</td><td>86,86%</td></tr><tr><td>LTC</td><td>100,000</td><td>5,12%</td><td>97,59%</td></tr><tr><td>LTC</td><td>500,000</td><td>0,24%</td><td>99,56%</td></tr></tbody></table>

A model with $\Delta(\theta) > 1$ performs worse than a simple polynomial fit under the prior.

Table 2 presents the values of $\Delta(\hat{\theta}_{norm})$ across the six datasets. As expected, discrepancies decrease with portfolio size. For annuities, the approximation performs reasonably well even at smaller scales. In contrast, for LTC, it yields substantial errors, except for the largest portfolio.

One explanation, supported by the Standardized Mortality Ratio (SMR) also provided in Table 2, is the positive correlation between observed event counts and their use as weights. This causes high crude rates to be overweighted, and low rates to be underweighted—introducing systematic overestimation of mortality rates. This bias is more severe in the LTC case where the observed deaths by data point is lower. In contrast, generalized WH smoothing preserves total event counts by construction, always yielding an SMR of exactly 100%. These results support adopting generalized WH smoothing in most practical settings. It retains the advantages of the original method while offering improved accuracy—even in small samples—and remains straightforward to implement.

# 5 How to select the smoothing parameters?

## 5.1 Impact of smoothing parameter choice

In the one-dimensional case, WH smoothing involves a single smoothing parameter $\lambda$; in two dimensions a pair $\lambda = (\lambda_x, \lambda_z)$. These parameters govern the trade-off between fidelity to the data and smoothness of the estimate, as defined in Equation 1.

Figure 1 illustrates this effect in a one-dimensional annuity dataset (100,000 policyholders, see Section 4.2), with three values of $\lambda$. The effective degrees of freedom (edf), computed as the trace of the hat matrix $H = (W + P_{\lambda})^{-1}W$, are shown for each curve. This quantity serves as a non-parametric analog of the number of free parameters in classical models and can take fractional values.
<caption>**Algorithm 1:** Iterative solution of generalized Whittaker-Henderson smoothing</caption>

inputs : d and e_c
outputs : ˆθ and Var(ˆθ)
parameters: λ, q, ε_l = 10⁻⁸
begin
    Construct the penalty matrix P_λ based on the difference matrices of order q.
    k ← 0
    θ₀ ← ln(d/e_c)
    l₀ ← -∞, cond_l ← true
    while cond_l do
        W_k ← Diag(exp(θ_k) ⊙ e_c)
        z_k ← θ_k + W_k⁻¹[d - exp(θ_k) ⊙ e_c]
        Ψ_{k+1} ← (W_k + P_λ)⁻¹
        θ_{k+1} ← Ψ_{k+1}W_kz_k
        l_{k+1} ← l_P(θ_{k+1}), cond_l ← l_{k+1} - l_k ≥ ε_l × sum(d)
        k ← k + 1
    ˆθ ← θ_k, Var(ˆθ) ← Ψ_k

For each of these portfolios, we first calculate the observed event counts **d** and the central exposure to risk **e**<sub>c</sub> based on discretized explanatory variables of interest: age *x* for the one-dimensional case and combinations (*x*, *z*) of age and long-term care duration for the two-dimensional case. We then apply Whittaker-Henderson smoothing to the vectors **y** = ln(**d**/**e**<sub>c</sub>) and **w** = **d** on one hand, or directly to the vectors **d** and **e**<sub>c</sub> within the framework of the generalization proposed in Section 4 on the other hand. The parameter λ used is determined for each example using a method that will be detailed in Section 5. Here, we simply specify that for each portfolio, the same parameter λ is used for both approaches, which corresponds to the same prior belief on the vector **θ** = ln **μ**.

The two approaches aiming to estimate the vector $\theta = \ln \mu$ can be directly compared using the parameters $\hat{\theta}_{\text{norm}}$ and $\hat{\theta}_{\text{ML}}$, estimated respectively by the normal approximation of the smoothing and by the maximization of the penalized likelihood. Since the first estimator can be seen as an approximate version of the second, we propose to use the error measure:

$$ \Delta(\theta) = \frac{\ell_P(\hat{\theta}_{\text{ML}}) - \ell_P(\theta)}{\ell_P(\hat{\theta}_{\text{ML}}) - \ell_P(\hat{\theta}_{\infty})} \quad (11) $$

where $\hat{\theta}_{\infty}$ corresponds to the parameter vector that maximizes the penalized likelihood with an infinite penalty. It will become apparent in Section 6.2 that this choice corresponds to the polynomial of degree $q-1$ that maximizes the likelihood. By definition, $\hat{\theta}_{\text{ML}} = \underset{\theta}{\operatorname{argmax}} \ell_P(\theta)$, and thus $\Delta(\theta) \ge 0$ for any vector $\theta$. Furthermore, $\Delta(\hat{\theta}_{\text{ML}}) = 0$ and $\Delta(\hat{\theta}_{\infty}) = 1$. A model
As shown, a low value $\lambda = 10^1$ yields an overfitted result that mirrors sampling noise, while a high value $\lambda = 10^7$ oversmooths and obscures the underlying trend. A mid-range value $\lambda = 10^4$ appears visually balanced. However, selecting a smoothing parameter by eye is unreliable: small-sample variability at the extremes of the age range can easily be mistaken for meaningful patterns.

<figure><img src="image_3.png" /><figcaption>Figure 1: WH smoothing on a synthetic annuity portfolio with 3 smoothing levels. Dots: crude rates; curves: smoothed estimates; shaded areas: credibility intervals. *edf*: effective degrees of freedom.</figcaption></figure>

The two-dimensional case further illustrates this difficulty. Figure 2 presents the smoothed transition rates from disability to death in an LTC portfolio (100,000 policyholders), using 9 combinations of ($\lambda_x, \lambda_z$). Choosing an appropriate parameter pair visually becomes nearly impossible, reinforcing the need for a data-driven statistical selection criterion.

## 5.2 Statistical criteria for parameter selection

Smoothing parameter selection typically relies on two classes of statistical criteria:

1. Prediction-based criteria, which aim to minimize prediction error, such as the Akaike Information Criterion (AIC) (Akaike 1973) and Generalized Cross-Validation (GCV) (Wahba 1980);

2. Likelihood-based criteria, which maximize the marginal likelihood—an approach introduced by Patterson and Thompson (1971) (under the name REML in the Gaussian case) and adapted to smoothing by Anderssen and Bloomfield (1974).

While prediction-based criteria have desirable asymptotic properties (Wahba 1985; Kauermann 2005), their convergence toward optimal smoothing parameters can be slow. In contrast,
<table><caption>Table 1: Impact of the normal approximation associated with Whittaker smoothing on a selection of 6 fictional insurance portfolios.</caption><thead><tr><td>Portfolio</td><td>Head count</td><td>Relative error on penalized deviance</td></tr></thead><tbody><tr><td>One-dimensional</td><td>20 000</td><td>7,52 %</td></tr><tr><td>One-dimensional</td><td>100 000</td><td>0,23 %</td></tr><tr><td>One-dimensional</td><td>500 000</td><td>0,09 %</td></tr><tr><td>Two-dimensional</td><td>1 000</td><td>1 784,54 %</td></tr><tr><td>Two-dimensional</td><td>5 000</td><td>122,87 %</td></tr><tr><td>Two-dimensional</td><td>25 000</td><td>6,95 %</td></tr></tbody></table>

for which $\Delta(\theta) > 1$ can be considered as having no practical interest, as it is less probable than a simple polynomial fit according to the prior belief used. The indicator $\Delta$ constructed in this way allows interpreting the quality of the used approximation within certain limits. For the purposes of this paper, a difference of less than 1% is considered negligible, while a difference greater than 10% is considered prohibitive. Differences between 1% and 10% will be interpreted more nuanced. This choice is obviously arbitrary.

The discrepancies obtained for each of the 6 portfolios studied are presented in Table 1. These discrepancies naturally decrease with the size of the portfolio, as the validity of the normal approximation increases. In the one-dimensional case, the discrepancy is modest for the two largest portfolios but significant for the smallest one. In the two-dimensional case, the discrepancy is significant for the largest portfolio and prohibitive for the two smallest ones. This is because using the observed death vector as the weight vector introduces a bias by overweighting observations with the highest crude rate and underweighting those with the lowest crude rate. In light of these results, it seems appropriate to prioritize the use of the generalized Whittaker-Henderson smoothing, especially since its implementation does not pose significant practical difficulties.

# 5 How to choose the smoothing parameter(s)?

This section is dedicated to the choice of the smoothing parameter $\lambda$ in the one-dimensional case or the pair $\lambda = (\lambda_x, \lambda_z)$ in the two-dimensional case. As shown in Figures 1 and 2, this choice strongly impacts the obtained results. Ideally, the selection of the smoothing parameter is based on the optimization of a statistical criterion, typically belonging to one of two major families.

On one hand, there are criteria based on the minimization of the model's prediction error, among which the Akaike information criterion (AIC, Akaike 1973) and the generalized cross-validation (GCV, Wahba 1980) are included. The Bayesian information criterion (BIC, Schwarz 1978) has a similar form to AIC, although its theoretical justification is very different, and thus it can be considered part of this group.
<figure><img src="image_4.png" /><figcaption>Figure 2: WH smoothing applied to disability-to-death transitions in an LTC portfolio, using 9 combinations of smoothing parameters. Contour lines and colours show the smoothed mortality surface by age and LTC duration.</figcaption></figure>

marginal likelihood criteria tend to perform more robustly in finite samples (Reiss and Todd Ogden 2009; Wood 2011).

To illustrate this, we apply AIC, GCV, and marginal likelihood to 100 replicates of the annuity portfolio with 100,000 policyholders (see Section 4.2). For each replicate, we select the optimal smoothing parameter and compute the corresponding effective degrees of freedom.

As shown on the left side of Figure 3, marginal likelihood produces stable and coherent degrees of freedom across replicates, whereas AIC and especially GCV often yield overly complex models. On the right, we plot the GCV and marginal likelihood profiles for a single replicate: marginal likelihood exhibits a well-defined maximum, while GCV presents two local minima. One aligns with the marginal likelihood optimum, but the global minimum corresponds to a model with ~35 degrees of freedom—an implausibly complex mortality curve.
On the other hand, there are criteria based on the maximization of a likelihood function known as the marginal likelihood. This type of criterion was introduced by Patterson H. D. (1971) in the Gaussian case, initially under the name of restricted likelihood (REML), and used by Anderssen and Bloomfield (1974) for the selection of smoothing parameters. Wahba (1985) and Kauermann (2005) show that criteria minimizing prediction error have the best asymptotic performance, but their convergence to the optimal smoothing parameters is slower. For finite sample sizes, criteria based on the maximization of a likelihood function are considered a more robust choice by many authors, such as Reiss and Todd Ogden (2009) or Wood (2011).

Figure 3 shows the values of GCV and marginal likelihood as a function of the smoothing parameter for the data used to produce Figure 1. When the marginal likelihood exhibits a clearly defined optimum, GCV has two optima, one of which coincides with the marginal likelihood’s optimum. The second optimum results in a model with nearly 45 degrees of freedom, which is not a plausible choice for the underlying mortality risk. For these reasons, we will favor the marginal likelihood as the selection criterion, especially since this choice naturally fits into the Bayesian framework introduced in Sections 2 and 4.

We will first discuss the selection in the case of the original Whittaker-Henderson smoothing and then extend the selection to the generalized case introduced in Section 4, introducing two competing approaches that will be subsequently compared.

<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison, in the context of one-dimensional Whittaker-Henderson smoothing, of the Generalized Cross-Validation (GCV) criterion and the marginal likelihood, as selection criteria for the smoothing parameter.</figcaption></figure>

## 5.1 Selection in the context of the original smoothing

Let us start from the notations and assumptions from Section 2, namely $y|\theta \sim N(\theta, W^{-1})$ and $\theta|\lambda \sim N(0, P_{\lambda}^{-1})$. In a purely Bayesian approach, it would be necessary to define a prior
These observations support the use of marginal likelihood over prediction-based criteria, especially in actuarial applications where robustness is key. Moreover, this choice aligns naturally with the Bayesian framework introduced in Sections 2 to 4.

We now detail its implementation—first for the original WH smoothing, then for the generalized setting—introducing three optimization strategies and three numerical algorithms and comparing their respective performances.

<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison of criteria for selecting the smoothing parameter in one-dimensional WH smoothing. Left: distribution of effective degrees of freedom under AIC, GCV, and marginal likelihood across 100 replicates. Right: GCV and marginal likelihood values for one replicate as functions of the smoothing parameter.</figcaption></figure>

## 5.3 Selection in the original smoothing

We consider again the normal framework from Section 2, where $y | \theta \sim N(\theta, \sigma^2 W^{-})$ and $\theta | \lambda \sim N(0, \sigma^2 P_{\lambda}^{-})$. In the empirical Bayes approach, the smoothing parameter $\lambda$ is estimated by maximizing the marginal likelihood:

$$ \mathcal{L}_{\text{norm}}^{m}(\lambda) = f(\mathbf{y} | \lambda) = \int f(\mathbf{y}, \theta | \lambda) d\theta = \int f(\mathbf{y} | \theta) f(\theta | \lambda) d\theta. $$

This is simply the maximum likelihood method applied to the smoothing parameter, treated as deterministic but unknown. A closed-form expression for this integral can be derived using standard Gaussian identities (see Section C of the appendices), yielding the marginal log-likelihood:

$$ \ell_{\text{norm}}^{m}(\lambda) = -\frac{1}{2} \left[ (\mathbf{y} - \hat{\theta}_{\lambda})^{T} W (\mathbf{y} - \hat{\theta}_{\lambda}) / \sigma^2 + \hat{\theta}_{\lambda}^{T} P_{\lambda} \hat{\theta}_{\lambda} / \sigma^2 + \ln|W + P_{\lambda}| - \ln|P_{\lambda}|_{+} + C \right]. $$

where $\hat{\theta}_{\lambda} = (W + P_{\lambda})^{-1} W \mathbf{y}$, and $C = -\ln|W|_{+} + (n_{*} - q) \ln(2\pi\sigma^2)$ is a constant independent of $\lambda$. This function is maximized numerically to obtain $\hat{\lambda}_{\text{norm}}$.
distribution on $\lambda$ and then estimate the posterior distribution of each parameter vector using methods such as Markov Chain Monte Carlo. The empirical Bayesian approach we adopt seeks to find the value of $\lambda$ that maximizes the marginal likelihood:

$$ \mathcal{L}_{\text{norm}}^{m}(\lambda) = f(\mathbf{y}|\lambda) = \int f(\mathbf{y}, \boldsymbol{\theta}|\lambda)d\boldsymbol{\theta} = \int f(\mathbf{y}|\boldsymbol{\theta})f(\boldsymbol{\theta}|\lambda)d\boldsymbol{\theta}. $$

This corresponds to the maximum likelihood method applied to the smoothing parameter. Let us explicitly rewrite the expressions of $f(\mathbf{y}|\boldsymbol{\theta})$ and $f(\boldsymbol{\theta}|\lambda)$ introduced in Section 2:

$$ f(\mathbf{y}|\boldsymbol{\theta}) = \sqrt{\frac{|W|_+}{(2\pi)^{n_*}}} \exp\left(-\frac{1}{2}(\mathbf{y}-\boldsymbol{\theta})^T W (\mathbf{y}-\boldsymbol{\theta})\right) $$

$$ f(\boldsymbol{\theta}|\lambda) = \sqrt{\frac{|P_{\lambda}|_+}{(2\pi)^{p-q}}} \exp\left(-\frac{1}{2}\boldsymbol{\theta}^T P_{\lambda} \boldsymbol{\theta}\right) $$

where $|A|_+$ denotes the product of the non-zero eigenvalues of $A$, $n_*$ is the number of non-zero diagonal elements of $W$, and $q$ is the number of zero eigenvalues of $P_\lambda$ ($q = q_x \times q_z$ in the two-dimensional case). Based on the Taylor expansion performed in Section 2, let us recall that:

$$ \ln f(\mathbf{y}, \boldsymbol{\theta}|\lambda) = \ln f(\mathbf{y}, \hat{\boldsymbol{\theta}}_\lambda|\lambda) + \frac{1}{2}(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_\lambda)^T (W + P_\lambda) (\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_\lambda) $$

which leads to:

$$ \begin{aligned} \mathcal{L}_{\text{norm}}^{m}(\lambda) &= \int \exp[\ln f(\mathbf{y}, \boldsymbol{\theta}|\lambda)] d\boldsymbol{\theta} \\ &= f(\mathbf{y}, \hat{\boldsymbol{\theta}}_\lambda|\lambda) \int \exp\left[-\frac{1}{2}(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_\lambda)^T (W + P_\lambda)(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_\lambda)\right] d\boldsymbol{\theta} \\ &= f(\mathbf{y}, \hat{\boldsymbol{\theta}}_\lambda|\lambda) \sqrt{\frac{(2\pi)^p}{|W + P_\lambda|}} \\ &= f_{\mathbf{y}}(\mathbf{y}|\hat{\boldsymbol{\theta}}_\lambda) f_{\boldsymbol{\theta}}(\hat{\boldsymbol{\theta}}_\lambda|\lambda) \sqrt{\frac{(2\pi)^p}{|W + P_\lambda|}} \\ &= \sqrt{\frac{|W|_+ |P_\lambda|_+}{(2\pi)^{n_*-q} |W + P_\lambda|}} \exp\left(-\frac{1}{2}\left[(\mathbf{y} - \hat{\boldsymbol{\theta}}_\lambda)^T W (\mathbf{y} - \hat{\boldsymbol{\theta}}_\lambda) + \hat{\boldsymbol{\theta}}_\lambda^T P_\lambda \hat{\boldsymbol{\theta}}_\lambda\right]\right). \end{aligned} $$

The associated log-likelihood can be expressed as follows:

$$ \begin{aligned} \ell_{\text{norm}}^{m}(\lambda) = & -\frac{1}{2} \left[ (\mathbf{y} - \hat{\boldsymbol{\theta}}_{\lambda})^{T} W (\mathbf{y} - \hat{\boldsymbol{\theta}}_{\lambda}) + \hat{\boldsymbol{\theta}}_{\lambda}^{T} P_{\lambda} \hat{\boldsymbol{\theta}}_{\lambda} \right. \\ & \left. \qquad - \ln|W|_{+} - \ln|P_{\lambda}|_{+} + \ln|W + P_{\lambda}| + (n_{*} - q) \ln(2\pi) \right]. \end{aligned} \qquad (12) $$
## 5.4 Selection in the generalized smoothing

The empirical Bayes approach introduced in the normal framework can be extended to the generalized smoothing framework developed in Section 4. While no closed-form expression exists for the marginal likelihood in this context, it can be approximated using a second-order Taylor expansion of the log-posterior density around its maximum $\hat{\theta}_{\lambda}$—similarly to what was done in the normal case. This yields the so-called Laplace Approximation of the Marginal Likelihood (LAML), defined as:

$$ \ell_{\text{LAML}}^{m}(\lambda) = \ell(\hat{\theta}_{\lambda}) - \frac{1}{2} \left[ \hat{\theta}_{\lambda}^{T} P_{\lambda} \hat{\theta}_{\lambda} + \ln|W_{\lambda} + P_{\lambda}| - \ln|P_{\lambda}|_{+} - q \ln(2\pi) \right] $$

where $W_{\lambda} = \text{Diag}(\exp(\hat{\theta}_{\lambda}) \odot e_{c})$ and $\ell(\hat{\theta}_{\lambda})$ is the log-likelihood evaluated at the penalized MLE. The detailed derivation of the Laplace approximation in this setting is provided in Section C of the appendices. This approximation plays a central role in the automatic selection of the smoothing parameter $\lambda$ in the generalized Whittaker-Henderson smoothing framework. As in the normal case, the marginal likelihood $\ell_{\text{LAML}}^{m}(\lambda)$ must be maximized numerically. However, a key distinction is that the penalized likelihood maximizer $\hat{\theta}_{\lambda}$ now depends on $\lambda$ and must be recomputed at each iteration via the PIRLS algorithm. This leads to a two-level optimization procedure:

* an inner loop estimating $\hat{\theta}_{\lambda}$ for fixed $\lambda$ using PIRLS;

* and an outer loop optimizing $\ell_{\text{LAML}}^{m}(\lambda)$ with respect to $\lambda$.

This outer iteration approach is the most principled method for smoothing parameter selection in this setting.

Alternative strategies have been proposed to reduce computational burden. The first one, known as performance-oriented iteration, was introduced by Gu (1992) and relies on the observation that, at each PIRLS step, the working response vector $z_k$ can be treated as approximately normal: $z_k | \theta \sim N(\theta, W_k^{-1})$. Assuming $W_k$ independent of $\lambda$, the marginal likelihood can be maximized within each PIRLS step using the normal approximation methodology of Section 5.3, with $y$ replaced by $z_k$ and $W$ by $W_k$. This effectively reverses the nesting structure, potentially saving computational time when updating $\lambda$ is less costly than recomputing a PIRLS step. A formal justification of the method is provided by Wood (2017, 149) which emphasizes that it does not actually require $z_k$ to have a normal distribution to be well-founded.

A third and even simpler strategy is the alternate iteration approach, used for instance by Wood et al. (2017). It consists in alternating updates of $\theta$ (via PIRLS) and $\lambda$ (via approximate marginal likelihood), without fully optimizing either at each step. This relies on the empirical observation that a coarse update of $\lambda$ may suffice, as the marginal likelihood surface changes between iterations.

Despite their efficiency, both performance-oriented and alternate iteration approaches lack formal convergence guarantees. Unlike outer iteration, they operate on different smoothing
Once the selection of $\lambda$ has been determined using the estimator $\hat{\lambda} = \underset{\lambda}{\operatorname{argmax}} \ell_{\text{norm}}^m(\lambda)$, the lack of an explicit solution to Equation 12 forces us to resort to numerical methods for its resolution. The Newton algorithm could once again be employed here and is a robust choice. This approach was notably adopted by Wood (2011). However, explicitly calculating the derivatives of the likelihood $\ell_{\text{norm}}^m$ is rather difficult from an operational perspective. Instead, we will use general heuristics such as those provided by Brent (1973) and John A. Nelder and Mead (1965), which are applicable to any sufficiently smooth function. These heuristics do not require derivative calculations and are implemented in the `optimize` and `optim` functions of the statistical programming language R.

Computing the marginal likelihood using Equation 12 requires the estimation of various intermediate quantities, with the most demanding being the estimation of $\hat{\theta}_{\lambda}$. This estimation is based on Equation 3 and ideally involves the inversion of the symmetric matrix $W + P_{\lambda}$ through a Cholesky decomposition. The calculation of $\ln|W + P_{\lambda}|$ is immediate for the triangular matrices resulting from this decomposition. Note that since the matrix $P_{\lambda}$ is a linear combination of the matrix $D_{n,q}^T D_{n,q}$ in the one-dimensional case, and the matrices $I_{n_z} \otimes D_{n_x,q_x}^T D_{n_x,q_x}$ and $D_{n_z,q_z}^T D_{n_z,q_z} \otimes I_{n_x}$ in the two-dimensional case, it is only necessary to form these matrices once. Furthermore, the calculation of $\ln|P_{\lambda}|_+$ can be efficiently performed by computing the eigendecomposition of the matrix $D_{n,q}^T D_{n,q}$ (or the matrices $D_{n_x,q_x}^T D_{n_x,q_x}$ and $D_{n_z,q_z}^T D_{n_z,q_z}$ in the two-dimensional case). This calculation then only requires multiplying these eigenvalues by the corresponding element of $\lambda$ and taking the logarithm. Finally, since the terms $\ln|W|_+$ and $(n_* - q) \ln(2\pi)$ do not depend on $\lambda$, they can be ignored altogether.

## 5.2 Selection in the generalized smoothing framework: outer iteration approach

It is possible to extend the previous approach to the framework of penalized likelihood introduced in Section 4. The Taylor expansion used in Section 5.1 can be applied in this context and yields:

$$ 
\begin{aligned}
\mathcal{L}_{\text{ML}}^{m}(\lambda) &= \int \exp[\ln f(\mathbf{d}, \mathbf{e}_{\mathbf{c}}, \boldsymbol{\theta}|\lambda)] d\boldsymbol{\theta} \\
&\approx f(\mathbf{d}, \mathbf{e}_{\mathbf{c}}, \hat{\boldsymbol{\theta}}_{\lambda}|\lambda) \int \exp\left(-\frac{1}{2}(\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_{\lambda})^T (W_{\hat{\boldsymbol{\theta}}_{\lambda}} + P_{\lambda}) (\boldsymbol{\theta} - \hat{\boldsymbol{\theta}}_{\lambda})\right) d\boldsymbol{\theta} \\
&\approx f(\mathbf{d}, \mathbf{e}_{\mathbf{c}}, \hat{\boldsymbol{\theta}}_{\lambda}|\lambda) \sqrt{\frac{(2\pi)^p}{|W_{\hat{\boldsymbol{\theta}}_{\lambda}} + P_{\lambda}|}} \\
&\approx l(\hat{\boldsymbol{\theta}}_{\lambda}) f_{\theta}(\hat{\boldsymbol{\theta}}_{\lambda}|\lambda) \sqrt{\frac{(2\pi)^p}{|W_{\hat{\boldsymbol{\theta}}_{\lambda}} + P_{\lambda}|}} \\
&\approx l(\hat{\boldsymbol{\theta}}_{\lambda}) \sqrt{\frac{(2\pi)^q |P_{\lambda}|_+}{|W_{\hat{\boldsymbol{\theta}}_{\lambda}} + P_{\lambda}|}} \exp\left(-\frac{1}{2} \hat{\boldsymbol{\theta}}_{\lambda}^T P_{\lambda} \hat{\boldsymbol{\theta}}_{\lambda}\right).
\end{aligned}
$$
parameters at each step, rendering penalized likelihood values non-comparable across iterations. Moreover, they do not track the value of $\ell_{LAML}^m(\lambda)$ during the optimization, making it harder to assess convergence or apply step-length controls.

Detailed algorithmic formulations of all three strategies in the generalized WH smoothing framework are provided in Section D of the appendices.

## 5.5 Algorithms for the maximization of the marginal likelihood

Several algorithms can be used to maximize the marginal likelihood or its Laplace approximation (LAML). It is generally preferable to apply these algorithms to the logarithm of the smoothing parameters, for three main reasons:

1. It ensures positivity of the smoothing parameters;

2. It simplifies the expressions of derivatives, when required;

3. It allows more uniform coverage of the range of interest (e.g., from $\lambda = 10^1$ to $10^7$, as in Figure 1, differences of comparable magnitude occur on a logarithmic scale).

### Derivative-free heuristics

A first, operationally simple option is to use general-purpose derivative-free optimization methods:

* Brent's method (Brent 1973) in the one-dimensional case;

* The Nelder-Mead simplex algorithm (Nelder and Mead 1965) in higher dimensions.

These are readily available in base R via the `optimize` and `optim` functions. They only require evaluating the marginal likelihood or LAML at each step, which is computationally inexpensive. However, they typically require more iterations to converge and cannot be combined with the alternate iteration approach, as they do not guarantee systematic improvement of the criterion at each step.

### Generalized Fellner-Schall method

A more specialized algorithm is the generalized Fellner-Schall method, based on ideas from Fellner (1986) and Schall (1991), and adapted for smoothing parameter selection in multidimensional generalized linear models by Rodriguez-Alvarez et al. (2015). It may be summarised by the update formula:

$$ \lambda_j^{\text{next}} = \frac{\text{tr}(P_{\lambda}^{-} P_j) - \text{tr}[(X^T W X + P_{\lambda})^{-1} P_j]}{\hat{\beta}_{\lambda}^T P_j \hat{\beta}_{\lambda}} \lambda_j^{\text{current}} \quad \text{for } j \in \{x, z\}. \quad (12) $$
This gives us the expression for the likelihood in this framework:

$$ \ell_{\text{ML}}^{m}(\lambda) \approx l(\hat{\theta}_{\lambda}) - \frac{1}{2} \left[ \hat{\theta}_{\lambda}^{T} P_{\lambda} \hat{\theta}_{\lambda} - \ln |P_{\lambda}|_{+} + \ln |W_{\hat{\theta}_{\lambda}} + P_{\lambda}| - q \ln(2\pi) \right]. \quad (13) $$

Unlike the normal case, the Taylor expansion does not allow for an exact calculation of the marginal likelihood but provides an approximation known as the Laplace approximation, whose validity depends on the number of available observations. Equation 13 does not have an explicit solution and requires numerical resolution. The calculation of the Laplace approximation of the marginal likelihood is done in a similar manner to Equation 12 from Section 5.1, with the only notable difference being that the vector $\hat{\theta}_{\lambda}$ must be estimated iteratively using Algorithm 1. This leads to Algorithm 2 for selecting the smoothing parameter which includes two nested iterative calculations: the calculation of $\hat{\theta}_{\lambda}$, for a fixed $\lambda$, using the Newton algorithm, and the search for the maximum of $\ell_{\text{ML}}^{m}$ using heuristics such as the Brent and Nelder-Mead algorithms in the one-dimensional and two-dimensional case, respectively. Since the selection of $\lambda$ corresponds to the outer iterative loop here, and to distinguish it from the approach introduced in the next section, we will refer to it as the outer iteration approach.

**Algorithm 2:** Parameter selection for generalized Whittaker-Henderson smoothing - outer iteration approach.

inputs : d and e_c
outputs : $\hat{\lambda}$, $\hat{\theta}_{\hat{\lambda}}$ and Var($\hat{\theta}_{\hat{\lambda}}$)
parameters: $q$, $\epsilon_{\ell} = 10^{-8}$, $\epsilon_{\text{ml}} = 10^{-8}$
begin
    k ← 0
    ml_0 ← -∞, cond_{ml} ← true
    while cond_{ml} do
        If $k = 0$, choose an arbitrary value $\lambda_0$ for the smoothing parameter(s); otherwise, choose the next value $\lambda_k$ using the selected heuristic
        Use Algorithm 1 to determine the vector $\theta_{\lambda_k}$ associated with the choice of $\lambda_k$ with an accuracy greater than $\epsilon_{\ell}$
        Calculate the marginal likelihood $\ell_{\text{ML}}^{m}(\lambda_k)$ associated with the choice of $\lambda_k$ using the intermediate quantities calculated during the estimation of $\theta_{\lambda_k}$
        $ml_{k+1} ← \ell_{\text{ML}}^{m}(\lambda_k)$, cond_{ml} ← $ml_{k+1} - ml_k ≥ \epsilon_{\text{ml}} \times \text{sum}(d)$
        $k ← k + 1$
    $\hat{\theta} ← \hat{\theta}_{\hat{\lambda}}$, Var($\hat{\theta}$) ← $\Psi_k$
where for WH smoothing $P_j = D_{n+,q}^T D_{n+,q}$ in the one-dimensional case and $P_x$ (resp. $P_z$) is the marginal penalization matrix $I_{n_z} \otimes D_{n_x,q_x}^T D_{n_x,q_x}$ (resp. $D_{n_z,q_z}^T D_{n_z,q_z} \otimes I_{n_x}$) in the two-dimensional case.This update can be interpreted more intuitively as:

$$ \hat{\beta}_{\lambda}^{T}(\lambda_{j}^{\text{next}} P_{j}) \hat{\beta}_{\lambda} = \text{tr}[P_{\lambda}^{-} \lambda_{j}^{\text{current}} P_{j} - (X^{T} W X + P_{\lambda})^{-1} \lambda_{j}^{\text{current}} P_{j}] $$

where the right-hand side corresponds to an effective degrees of freedom associated with $\lambda_j^{\text{current}} P_j$, and the left-hand side to a squared error, normalized by the updated penalty precision. This makes $\lambda_j^{\text{next}}$ resemble a REML-based estimator for the inverse variance. More details may be found in Rodríguez-Álvarez et al. (2019). This method:

* May be combined with any of the three iteration nesting schemes (outer, performance, alternate);

* Does not require explicit derivative computations;

* Converges toward an approximate maximum of LAML in the generalized case, since it ignores the dependence of $W$ on $\lambda$;

* Tends to take longer steps than EM-like algorithms (Dempster, Laird, and Rubin 1977), but shorter than Newton updates (see Wood and Fasiolo 2017, which also provides a thorough justification for the method).

## Newton algorithm

A third option is the Newton method, which involves computing both the first and second derivatives of the marginal likelihood (or LAML) with respect to $\ln \lambda$. Full derivations are provided in Wood (2011), which covers a more general case. The method applies in both the normal and generalized cases, but in the latter, derivative expressions are more complex due to the dependence of $W$ on $\lambda$. The Newton algorithm is fast and precise and applicable to all three nesting strategies. The downside is the operational complexity associated with this method, especially in the generalized case.

## 5.6 Performance comparison

Sections 5.4 and 5.5 introduced eight combinations of nesting strategies and optimization algorithms applicable to the generalized WH smoothing. We now assess the potential convergence issues and approximation errors associated with each of them.

This analysis is based on 100 replicates of the simulated annuity and LTC portfolios with 100,000 policyholders, as described in Section 4.2. For each replicate and each method combination, we compute the LAML at the selected smoothing parameters, and compare this value to the (approximate) optimal value obtained across all combinations, denoted $\hat{\lambda}_{\text{opt}}$.
## 5.3 Selection in the generalized smoothing framework: performance iteration approach

An alternative to the approach developed in Section 5.2 is to start from Algorithm 1 and notice that at each step of the algorithm, the solved equations coincide with those obtained by assuming that $z_k|\theta \sim \mathcal{N}(\theta, W_k^{-1})$ and $\theta \sim \mathcal{N}(0, P_\lambda^{-1})$. Therefore, it is possible to rely on the methodology proposed in Section 5.1 and numerically estimate, at the beginning of each iteration of Algorithm 1, the smoothing parameter that maximizes the marginal likelihood $\ell_{\text{norm}}^m$ in which $y$ would be replaced by $z_k$ and $W$ by $W_k$. This is what Algorithm 3 proposes. A justification for this approach is given by (Wood 2006, p149). From a practical point of view, this approach conceptually reverses the nesting order of the two types of iterations and can, in some cases, provide a considerable time gain. It was introduced by Gu (1992) as the “performance (oriented) iteration” approach. However, unlike the outer iteration approach, the convergence of the performance iteration approach cannot be guaranteed. Indeed, the penalized likelihood calculated in Algorithm 3 is not directly comparable between iterations because it is based on a different smoothing parameter and, therefore, a different prior. There are situations in which this algorithm may not converge. However, this difficulty seems to be more of a theoretical nature than a practical one for the optimization problem addressed in this paper.

<caption>**Algorithm 3:** Parameter selection for generalized Whittaker-Henderson smoothing - performance iteration approach</caption>

inputs : $d$ and $e_c$
outputs : $\hat{\lambda}$, $\hat{\theta}_{\hat{\lambda}}$ and $\text{Var}(\hat{\theta}_{\hat{\lambda}})$
parameters: $q$, $\epsilon_\ell = 10^{-8}$, $\epsilon_{\text{ml}} = 10^{-8}$
begin
    $k \leftarrow 0$
    $\theta_0 \leftarrow \ln(d/e_c)$
    $l_0 \leftarrow -\infty$, $\text{cond}_\ell \leftarrow \text{true}$
    while $\text{cond}_\ell$ do
        $W_k \leftarrow \text{Diag}(\exp(\theta_k) \odot e_c)$
        $z_k \leftarrow \theta_k + W_k^{-1}[d - \exp(\theta_k) \odot e_c]$
        Estimate the parameter $\lambda_k$ maximizing the regularized marginal likelihood $\ell_{\text{norm}}^m$ associated with the observation vector $zk$ and the weight matrix $W_k$, using the selected heuristic, with an accuracy greater than $\epsilon_{\text{ml}} \times \text{sum}(d)$
        $\Psi_{k+1} \leftarrow (W_k + P_{\lambda_k})^{-1}$
        $\theta_{k+1} \leftarrow \Psi_{k+1} W_k z_k$
        $l_{k+1} \leftarrow \ell_P(\theta_{k+1}|\lambda_k)$, $\text{cond}_\ell \leftarrow l_{k+1} - l_k \ge \epsilon_\ell \times \text{sum}(d)$
        $k \leftarrow k + 1$
    $\hat{\lambda} \leftarrow \lambda_k$; Use Algorithm 1 to obtain the values of $\hat{\theta}_{\hat{\lambda}}$ and $\text{Var}(\hat{\theta}_{\hat{\lambda}})$ corresponding to the parameter $\hat{\lambda}$
To quantify the discrepancy, we define the relative error:

$$ \Delta(\lambda) = \frac{\ell_{\text{LAML}}^{m}(\hat{\lambda}_{\text{opt}}) - \ell_{\text{LAML}}^{m}(\lambda)}{\ell_{\text{LAML}}^{m}(\hat{\lambda}_{\text{opt}}) - \ell_{\text{LAML}}^{m}(\infty)} \quad (13) $$

where $\ell_{\text{LAML}}^{m}(\infty)$ corresponds to the LAML value when using an infinite smoothing penalty, i.e., the overly smooth baseline. By construction, $\Delta(\lambda) \ge 0$ for all tested methods, with $\Delta(\hat{\lambda}_{\text{opt}}) = 0$ and $\Delta(\infty) = 1$. In the two-dimensional setting, we also compare average computation time for each method across replicates.

Results are summarised in Figure 4. The top panel displays the relative error $\Delta(\lambda)$ (capped below $10^{-10}$ for readability). In the outer iteration framework:

* The Newton method consistently achieves relative errors below $10^{-10}$;

* Brent and Nelder-Mead heuristics yield slightly higher errors but remain below $10^{-7}$;

* The generalized Fellner-Schall method produces higher errors, but still below $10^{-5}$ and negligible in practice.

In the performance and alternate iteration frameworks, all methods yield similar errors, consistently below $10^{-5}$, with no convergence issues observed in any replicate. These findings suggest that method selection can be guided by practical considerations such as speed and implementation ease.

The bottom panel of Figure 4 compares computation times (relative to the Nelder-Mead + outer iteration baseline):

* In the outer iteration framework, the Newton method is the fastest, followed by the Fellner-Schall approach;

* All outer iteration variants are faster than their performance or alternate counterparts.

This is unsurprising, as PIRLS steps are particularly lightweight in WH smoothing (where the model matrix is the identity). However, alternate strategies may remain useful for more general cases like those described in Section 6.4.

For reference, the average time required for a single iteration using Nelder-Mead in the 2D outer iteration case is approximately 1.68 seconds (versus 5 milliseconds in the 1D case).
## 5.4 Comparison of outer iteration and performance iteration approaches

In Sections 5.2 and 5.3, we introduced two alternatives for the selection of parameters in generalized Whittaker-Henderson smoothing. This section aims to study the consequences of choosing one approach over the other in terms of result accuracy and computational time. While the objective of Section 4.2 was to highlight a bias in the asymptotic approximation associated with the original Whittaker-Henderson smoothing for small portfolios, here we aim to measure these impacts more precisely. We limit our analysis to intermediate-sized portfolios as presented in Section 4.2 (i.e., 100,000 rows for the one-dimensional case and 5,000 rows for the two-dimensional case), but we use 100 replicates of each portfolio, generated from the same mortality and censoring laws, to increase the robustness of our analysis.

The outer iteration and performance iteration approaches can be directly compared based on the selected parameter $\lambda$ and using the marginal likelihood $\ell_{ML}^m(\lambda)$. Similar to Section 4.2, we define the error criterion as:

$$ \Delta(\lambda) = \frac{\ell_{ML}^m(\hat{\lambda}_{outer}) - \ell_{ML}^m(\lambda)}{\ell_{ML}^m(\hat{\lambda}_{outer}) - \ell_{ML}^m(\infty)} \quad (14) $$

where $\ell_{ML}^m(\infty)$ corresponds to the approximation of the marginal likelihood associated with the choice of an infinite smoothing parameter. Again, the defined $\Delta$ satisfies the properties $\Delta(\lambda_{outer}) = 0$ and $\Delta(\infty) = 1$. This criterion is interpreted in the same way as the one given by Equation 11.

<figure><img src="image_6.png" /><figcaption>Figure 4: Distribution of the difference between the two presented approaches for selecting the smoothing parameter, for 100 simulated portfolios, in the one-dimensional and two-dimensional case.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Comparison of the 8 nesting strategy and algorithm combinations in the 1D and 2D simulated cases. Top: relative error on the LAML (log scale). Bottom: improvement in average computation time compared to the Nelder-Mead + outer iteration reference.</figcaption></figure>

# 6 How to improve smoothing computational efficiency?

## 6.1 Motivation

Whittaker-Henderson (WH) smoothing is a full-rank method, meaning that it includes as many parameters as there are observation points. This feature ensures a high degree of flexibility, allowing the estimator to closely track the input signal when sufficient data is
Figure 4 represents the empirical distribution of the difference between the outer iteration and performance iteration approaches, and Figure 5 represents the distribution of the associated computation time for each approach. The use of the performance iteration approach results in a difference below 1% for the 100 simulated datasets, both in the one-dimensional and two-dimensional case. The associated computation time is reduced by a factor of 2 in the one-dimensional case and nearly 4 in the two-dimensional case. Considering the significant computation times in the latter case, the use of the performance iteration approach presents a real operational advantage.

<figure><img src="image_7.png" /><figcaption>Figure 5: Distribution of the computation time associated with the two presented approaches for selecting the smoothing parameter, for 100 simulated portfolios, in the one-dimensional and two-dimensional case.</figcaption></figure>

# 6 How to improve smoothing performances with a large number of data points?

## 6.1 Motivation

Whittaker-Henderson smoothing is considered to be full-rank smoothing as it contains as many parameters as there are observation points. This characteristic allows it to faithfully reproduce any input signal, provided that a sufficient number of observations are available. More formally, the estimator associated with smoothing is asymptotically unbiased since:

$$ \mathbb{E}(\hat{\mathbf{y}}) = (W + P_{\lambda})^{-1} W \mathbb{E}(\mathbf{y}) \underset{m \to \infty}{\to} \mathbb{E}(\mathbf{y}) $$

where m represents the number of observed individuals and is involved in the matrix W. The downside is that it may become impractical in presence of many observations points n. Indeed,
available. Formally, WH smoothing is asymptotically unbiased since:

$$ \mathbb{E}(\hat{\mathbf{y}}) = (W + P_{\lambda})^{-1} W \mathbb{E}(\mathbf{y}) \xrightarrow{m \to \infty} \mathbb{E}(\mathbf{y}), $$

where $m$ denotes the number of observed individuals, which influences the matrix $W$.

However, this flexibility comes at a computational cost. Some key operations, such as (implicit) matrix inversions, scale cubically with the number of parameters. As a result, WH smoothing may become impractical with large number of combinations or when applied repeatedly (e.g. in simulations or bootstraps).

In one-dimensional settings, such as age-only models with annual discretization, the number of points rarely exceeds 100, and computation time is negligible. In contrast, two-dimensional use cases—common in insurance—can lead to substantially larger datasets:

* Disability tables in France must cover entry ages from 18 to 61 and exit ages up to 62, resulting in $(62 - 18) \times (62 - 18 + 1)/2 = 990$ combinations.

* Transition tables from short-term incapacity to disability involve entry ages from 18 to 67 and monthly durations from 0 to 36 months, yielding $(67 - 18) \times 36 = 1,764$ combinations.

* Long-term care (LTC) models require coverage over ages 50 to 100 and durations from 0 to 20 years, totalling $(100 - 50) \times (20 - 0) = 1,000$ combinations (in practice, this number may be lower due to data sparsity).

In such settings, computing WH smoothing—especially when paired with smoothing parameter selection—can take several minutes per application, limiting usability in iterative contexts.

To address this limitation, we now analyse the computational complexity of the main steps in WH smoothing and smoothing parameter selection then introduce two complementary strategies to reduce computation time:

* A structural optimization that exploits the specific form of WH penalization matrices;

* A reduced-rank approximation that lowers the number of parameters while minimizing bias compared to the full-rank estimator.

Finally, we benchmark these strategies in terms of runtime and accuracy using 100 replicates of the mid-size annuity and LTC portfolios from Section 4.2. The structural optimization is compared to the original WH method, while the reduced-rank approximation is evaluated against the original method, the structural optimization, and a reference P-spline smoothing approach.
the algorithms presented in Section 5 require the inversion of the $n \times n$ matrix $W + P_\lambda$, an operation with a time complexity of $O(n^3)$, which needs to be repeated at each iteration of the algorithm, regardless of the chosen method.

For biometric risks that depend solely on age, discretized on an annual basis, the number of observations rarely exceeds 100, and computation time is not a significant concern. However, in the two-dimensional case, the number of observations can take much larger values in several practical cases, including:

* For the disability risk, it is necessary to construct disability survival tables for entry ages ranging from 18 to 62 and exit ages ranging from the entry age to 62. This represents $[(62 - 18) \times (62 - 18 + 1)/2 = 990$ observations.

* Also for the disability risk, it is necessary to construct transition tables from incapacity to disability for entry ages ranging from 18 to 67 and monthly incapacity durations ranging from 0 to 36 months. This represents $(67 - 18) \times (36 - 0) = 1,764$ observations.

* For the long-term care risk, as the coverage is lifelong, it is necessary to construct tables for all ages and long-term care durations encountered in the data. Assuming observation ages ranging from 50 to 110 and durations ranging from 0 to 20 years, this already represents $(110 - 50) \times (20 - 0) = 1,200$ observations.

The practical applications described above thus require computation times on the order of several minutes for each smoothing application. When smoothing is intended to be used repeatedly, for example, in simulation contexts, reducing the computation time becomes an operational challenge.

## 6.2 Smoothing and eigendecomposition

The eigendecomposition of the penalization matrix is key to a better understanding of smoothing and allows for an approximation of the smoothing problem with $p < n$ parameters. Let us consider the one-dimensional case and write the decomposition for the symmetric matrix $D_{n,q}^T D_{n,q}$. It takes the form $D_{n,q}^T D_{n,q} = U \Sigma U^T$, where $\Sigma$ is a diagonal matrix containing the eigenvalues of $D_{n,q}^T D_{n,q}$ and $U$ is an orthogonal matrix such that $U^T U = U U^T = I_n$. Let us perform the reparameterization $\beta = U^T \theta \Leftrightarrow \theta = U\beta$. The smoothness criterion becomes:

$$ \theta^T P_\lambda \theta = (U\beta)^T P_\lambda (U\beta) = \lambda \beta^T U^T U \Sigma U^T U \beta = \lambda \beta^T \Sigma \beta. \quad (15) $$

and Equation 2 takes the form:

$$ \hat{y} = U\hat{\beta} \quad \text{where} \quad \hat{\beta} = \underset{\beta}{\operatorname{argmin}} \left\{ (\mathbf{y} - U\beta)^T W (\mathbf{y} - U\beta) + \lambda \beta^T \Sigma \beta \right\}. \quad (16) $$
## 6.2 Practical computation for penalized smoothers

Whittaker-Henderson (WH) smoothing belongs to a broader family of penalized smoothing methods that produce estimates of the form:

$$ \hat{y}_{\lambda} = X\hat{\beta}_{\lambda} \quad \text{where } \beta_{\lambda} \text{ solves } (X^{T}WX + P_{\lambda})\hat{\beta}_{\lambda} = X^{T}Wy. $$

Here, $X$ and $P_{\lambda}$ denote the model and penalization matrices of size $n \times p$ and $p \times p$ respectively, and $W$ is a diagonal matrix of positive weights of size $n \times n$.

### Computational steps

The computation of $\hat{y}_{\lambda}$ for a given $\lambda$ typically involves the following steps:

1. Absorb the weights in the model matrix and observation vector, forming $W^{1/2}X$ and $W^{1/2}y$, which requires $O(n^2)$ and $O(n)$ operations respectively (multiplying each row of $X$ and each element of $y$ by the corresponding element of $w$).

2. Form the matrix $P_{\lambda}$. The cost of this operation is typically $O(p^2)$ in the general case.

3. Form the matrix $X^{T}WX$ and the vector $X^{T}Wy$, which requires up to $O(np^2)$ and $O(np)$ operations respectively.

4. Add together $X^{T}WX$ and $P_{\lambda}$ which requires $O(p^2)$ operations in the general case.

5. Compute the Cholesky decomposition $X^{T}WX + P_{\lambda} = R^{T}R$ at a cost of $O(p^3)$.

6. Obtain $\hat{\beta}_{\lambda}$ by forward-backward substitution, first solving $R^{T}u = X^{T}Wy$ then $R\hat{\beta}_{\lambda} = u$ with an associated cost of $O(p^2)$ for each system.

7. Compute $\hat{y}_{\lambda} = X\hat{\beta}_{\lambda}$ at a cost of $O(np)$.

As an alternative to Cholesky, QR decomposition may be used for greater numerical stability (see Golub and Van Loan 2013). It applies to the weighted design matrix stacked with a matrix $B$ such that $B^{T}B = P_{\lambda}$.

### Simplifications for WH smoothing

In WH smoothing, $X = I_n$, which simplifies computations:

* Step 7 is unnecessary, as well as the first part of step 3.

* $X^{T}Wy = Wy$ (step 3) is computed in $O(n)$ by multiplying $w$ and $y$.

* $X^{T}WX + P_{\lambda} = W + P_{\lambda}$ (step 4) is also computed in $O(n)$ by adding the vector $w$ to the leading diagonal of $P_{\lambda}$.
In the original formulation of smoothing, the parameters directly corresponded to the smoothed values. With this new parameterization, $\beta$ can be interpreted as a vector of coordinates in the basis of eigenvectors of $P_\lambda$, providing a decomposition of the signal into components that are more or less smooth according to the employed penalization. Figure 6 represents 8 of the eigenvectors associated with $q = 2$ for a basis of size $n = 74$, which corresponds to the number of observations in the mortality portfolio used, notably in Figure 1. The eigenvalues associated with the first $q$ eigenvalues are zero. This result can be easily explained by noting that $D_{n,q}$ is a matrix of dimensions $(n-q) \times n$ and rank $n-q$.

<figure><img src="image_8.png" /><figcaption>Figure 6: Graphical representation of a subset of columns from a basis of eigenvectors associated with $n = 74$ and $q = 2$ in the one-dimensional case. The 8 vectors in the subset are displayed in order of increasing penalization.</figcaption></figure>

By using the fact that $U^{-1} = U^T$ and making the connection with Equation 3, we obtain the explicit solution:

$$ \hat{y} = U(U^T W U + S_\lambda)^{-1} U^T W y \quad \text{where} \quad S_\lambda = \lambda \Sigma \qquad (17) $$

In order to interpret Equation 17, let us consider the special case where all weights are equal to 1, and we have:
## Generalized WH smoothing with outer iteration

When using the outer iteration approach (see Section 5), each candidate $\lambda$ requires a full PIRLS cycle to estimate $\hat{\theta}_{\lambda}$, with new working vector $\hat{z}_{\lambda}^k$ and weight matrix $W_{\lambda}^k$. Steps 1–6 above are repeated until convergence of the PIRLS algorithm, which may be assessed by monitoring the changes in penalized deviance. The deviance may be computed at a $O(n)$ cost. For penalization based on differences matrices, computation of $\hat{\beta}_{\lambda}^T P_{\lambda} \hat{\beta}_{\lambda}$ should be based on the expression of $R_{\lambda,q}$ provided in Section 1.1 for an associated cost of $O(qp)$. In addition, PIRLS iterations for each new $\lambda$ can be initialized using the previous estimate of $\hat{y}_{\lambda}$ for faster convergence.

## LAML computation

Once the deviance is known, computing the marginal likelihood/LAML also requires:

* $\ln|X^TWX + P_{\lambda}|$, which may be computed at a cost of $O(p)$ from the leading diagonal of the Cholesky/QR factor $R$ computed at step 5 in the derivation of $\hat{y}_{\lambda}$.

* $\ln|P_{\lambda}|_+$, which may be obtained from the eigenvalues of the penalization matrix: Section 6.4 shows that in the two-dimensional case, it can be computed via eigendecomposition of $D_{p_x,q_x}^T D_{p_x,q_x}$ and $D_{p_z,q_z}^T D_{p_z,q_z}$, performed only once, at a $O(p_x^3+p_z^3)$ cost. Computation of $\ln|P_{\lambda}|_+$ then only requires scaling the eigenvalues for a cost of $O(p)$.

## Algorithm-specific computations

Brent and Nelder-Mead require only marginal likelihood/LAML evaluations.

The generalized Fellner-Schall algorithm relies on the update formula of Equation 12:

$$ \lambda_j^{\text{next}} = \frac{\text{tr}(P_{\lambda}^{-} P_j) - \text{tr}[(X^T W X + P_{\lambda})^{-1} P_j]}{\hat{\beta}_{\lambda}^{T} P_j \hat{\beta}_{\lambda}} \quad \text{for } j \in \{x, z\}. $$

Evaluation of $\text{tr}(P_{\lambda}^{-} P_j)$ does not require any matrix product. In the one-dimensional case it is simply $(p-q)/\lambda$ while in the two-dimensional case it may be obtained directly at a $O(p)$ cost using the eigenvalues of the aforementioned penalization matrices $P_j$. Evaluation of $\text{tr}[(X^T W X + P_{\lambda})^{-1} P_j]$ may use the identity $\text{tr}(AB) = \sum_{i,j} A_{ij}B_{ji}$ and therefore be computed at an $O(p^2)$ cost if the matrix $(X^T W X + P_{\lambda})^{-1}$ and $P_j$ are available. Computation of $V = (X^T W X + P_{\lambda})^{-1}$ is done by first solving for the inverse $K = R^{-1}$ of the Cholesky/QR factor and then forming $V$ as $KK^T$. Both operations have a $O(p^3)$ cost.

Newton method also requires computation of $V$, as well as several matrix products involving the penalization matrix $P_j$. For example, the second derivatives of marginal likelihood require $\text{tr}[VP_jVP_k]$ terms and the second derivatives of marginal likelihood require $\text{tr}[V(X(\partial W/\partial \rho_j)X + P_j)V(X(\partial W/\partial \rho_k)X + P_k)]$ terms where $\rho_j = \ln(\lambda_j)$, $j = k = x$ in the one-dimensional case
$$ \hat{\mathbf{y}} = U(U^T U + \lambda \Sigma)^{-1} U^T \mathbf{y} = U(I_n + \lambda \Sigma)^{-1} U^T \mathbf{y} $$

The transformation from $\mathbf{y}$ to $\hat{\mathbf{y}}$ can then be seen as a three-step process, reading the equation from right to left:

1. Decomposition of the signal $\mathbf{y}$ in the basis of eigenvectors through the left multiplication by $U^T$.

2. Attenuation of the signal components based on the eigenvalues associated with these components. If we denote $s = \text{diag}(\Sigma)$, then $(I_n + \lambda \Sigma)^{-1} = \text{Diag}[1/(1 + \lambda s)]$. After the left multiplication by $(I_n + \lambda \Sigma)^{-1}$, each component is divided by a coefficient $1 + \lambda s \ge 1$. This attenuation coefficient increases linearly with $\lambda$, but at different rates for each eigenvalue.

3. Reconstruction of the attenuated signal in the canonical basis through the left multiplication by $U$.

In the presence of non-unit weights, things are not as straightforward since $U^T W U$ is no longer a diagonal matrix. However, it is still possible to visualize the effect of smoothing by looking at the matrix $F = (U^T W U + S_\lambda)^{-1} U^T W U$. In fact:

$$ U^T \hat{\mathbf{y}} = U^T U \hat{\mathbf{\theta}} = \hat{\mathbf{\theta}} = (U^T W U + S_\lambda)^{-1} U^T W \mathbf{y} = (U^T W U + S_\lambda)^{-1} U^T W U U^T \mathbf{y} = F U^T \mathbf{y} $$

Since the vectors $U^T \mathbf{y}$ and $\hat{\mathbf{\theta}} = U^T \hat{\mathbf{y}}$ represent the coordinates of the signals $\mathbf{y}$ and $\hat{\mathbf{y}}$ in the basis of eigenvectors of $D_{n,q}^T D_{n,q}$, $F$ corresponds to a coordinate transformation matrix playing a similar role for the parameters as the hat matrix $H = U(U^T W U + S_\lambda)^{-1} U^T W$ does for the observations. The diagonal values of $F$ can be interpreted as the effective degrees of freedom associated with each eigenvector after smoothing. It can be verified that:

$$ \text{tr}(F) = \text{tr}[(U^T W U + S_\lambda)^{-1} U^T W U] = \text{tr}[U(U^T W U + S_\lambda)^{-1} U^T W] = \text{tr}(H) $$

which means that the sum of the effective degrees of freedom remains the same whether it is counted per observation or per parameter.

Figure 7 represents the effective degrees of freedom per parameter in the previous illustration of smoothing. The first $q$ eigenvectors are never penalized, so their effective degrees of freedom are always equal to 1, regardless of the smoothing parameter used. The other eigenvectors have strictly decreasing effective degrees of freedom with $\lambda$. These degrees of freedom are generally decreasing with increasing eigenvalues of $D_{n,q}^T D_{n,q}$, although in the presence of non-unit weights and for small values of $\lambda$, this may not always be the case.
and {$j, k \in \{x, z\}$} in the two-dimensional case. The identity $\text{tr}(AB) = \sum_{i,j} A_{ij}B_{ji}$ can also be used in this case but matrix products $VP_j$ or $V[X(\partial W/\partial \rho_j)X + P_j]$ still need to be explicitly computed, for a respective cost of $O(p^3)$ and $O(np^2)$ each.

These additional computations make Newton updates more expensive than generalized Fellner-Schall updates, but they generally yield faster convergence and higher precision (see Section 5.6).

## 6.3 Banded optimization for WH smoothing

We now consider how to exploit the banded structure of the penalization matrix in Whittaker-Henderson (WH) smoothing. This structure enables significant computational gains, especially when dealing with large number of observations. Throughout this section, we assume $X = I_n$ and $p = n$, which holds for both the original and generalized WH smoothing.

### One-Dimensional Case

In one dimension, the penalization matrix takes the form: $P_\lambda = \lambda D_{n,q}^T D_{n,q}$. This matrix is symmetric and banded with bandwidth $q$. As a consequence:

1. Compact storage: $P_\lambda$ can be stored in a compact form with dimensions $n \times (q + 1)$, and updated for new $\lambda$ at a cost of $O(qn)$. The matrix $W + P_\lambda$ shares this structure.

2. Efficient Cholesky decomposition: the Cholesky factor $R$ of $W + P_\lambda$ can be computed in $O(q^2n)$ instead of $O(n^3)$, and $R$ is also banded with the same bandwidth.

3. Efficient back-substitution: computing $\hat{y}_\lambda = \hat{\beta}_\lambda$ using $R$ is now $O(qn)$ instead of $O(n^2)$.

4. Efficient inversion of $R$: the inverse $K = R^{-1}$ costs $O(qn^2)$, an improvement over the $O(n^3)$ cost for dense matrices.

However, $K$ is a dense triangular matrix, meaning the computation of $V = KK^T$ remains a $O(n^3)$ operation. Fortunately, the generalized Fellner-Schall algorithm only requires the diagonal of $V$, which can be obtained from $K$ in $O(n^2)$, since:

$$ \lambda_j[\text{tr}(P_\lambda^{-} P_j) - \text{tr}(VP_j)] = (n - q) - (n - \text{tr}[VW]) = \text{diag}(V)^T w - q. $$

Furthermore, the Newton algorithm benefits as well: the trace terms involving $VP_\lambda$ or $V(\partial W/\partial \ln \lambda + P_\lambda)$ are based on banded matrices, making those products computable in $O(qn^2)$ instead of $O(n^3)$.
<figure><img src="image_9.png" /><figcaption>Figure 7: Distribution of residual degrees of freedom per eigenvector after applying the Whittaker-Henderson smoothing, in the one-dimensional case, for different values of the smoothing parameter.</figcaption></figure>

## 6.3 Extension to the two-dimensional case

In the two-dimensional case, we have $P_{\lambda} = \lambda_{x}I_{n_{x}} \otimes D_{n_{x},q_{x}}^{T}D_{n_{x},q_{x}} + \lambda_{z}D_{n_{z},q_{z}}^{T}D_{n_{z},q_{z}} \otimes I_{n_{x}}$. Similar to the one-dimensional case, we can perform the eigendecomposition of the matrices $D_{n_{x},q_{x}}^{T}D_{n_{x},q_{x}}$ and $D_{n_{z},q_{z}}^{T}D_{n_{z},q_{z}}$, yielding $D_{n_{x},q_{x}}^{T}D_{n_{x},q_{x}} = U_{x}\Sigma_{x}U_{x}^{T}$ and $D_{n_{z},q_{z}}^{T}D_{n_{z},q_{z}} = U_{z}\Sigma_{z}U_{z}^{T}$. Let us define $U = U_{z} \otimes U_{x}$ and perform the reparameterization $\beta = U^{T}\theta \Leftrightarrow \theta = U\beta$. By leveraging the properties of the Kronecker product, we can rewrite the smoothness criterion in a simplified form:

$$ 
\begin{aligned}
\theta^T P_{\lambda} \theta &= (U\beta)^T P_{\lambda}(U\beta) \\
&= (U\beta)^T (\lambda_x I_{n_z} \otimes U_x \Sigma_x U_x^T + \lambda_z U_z \Sigma_z U_z^T \otimes I_{n_x}) (U\beta) \\
&= \beta^T (U_z \otimes U_x)^T (\lambda_x I_{n_z} \otimes U_x \Sigma_x U_x^T + \lambda_z U_z \Sigma_z U_z^T \otimes I_{n_x}) (U_z \otimes U_x)\beta \\
&= \beta^T (U_z^T \otimes U_x^T) (\lambda_x I_{n_z} \otimes U_x \Sigma_x U_x^T + \lambda_z U_z \Sigma_z U_z^T \otimes I_{n_x}) (U_z \otimes U_x)\beta \\
&= \beta^T (\lambda_x U_z^T U_z \otimes U_x^T U_x \Sigma_x U_x^T U_x + \lambda_z U_z^T U_z \Sigma_z U_z^T U_z \otimes U_x^T U_x)\beta \\
&= \beta^T (\lambda_x I_{n_z} \otimes \Sigma_x + \lambda_z \Sigma_z \otimes I_{n_x})\beta.
\end{aligned}
$$

This leads to an alternative formulation of the optimization problem:

$$ \hat{y} = U\hat{\beta} \quad \text{where} \quad \hat{\beta} = \underset{\beta}{\operatorname{argmin}} \left\{ (y - U\beta)^T W (y - U\beta) + \lambda\beta^T (\lambda_x I_{n_z} \otimes \Sigma_x + \lambda_z \Sigma_z \otimes I_{n_x}) \beta \right\}. \qquad (18) $$
## Two-dimensional Case

In two dimensions, the penalization matrix is $P_{\lambda} = \lambda_{x}P_{x} + \lambda_{z}P_{z}$ where:

* $P_x = I_{n_z} \otimes D_{n_x, q_x}^T D_{n_x, q_x}$

* $P_z = D_{n_z, q_z}^T D_{n_z, q_z} \otimes I_{n_z}$.

This structure has the following key properties:

* Both matrices are made of $n_z \times n_z$ square blocks of dimensions $n_x \times n_x$ each.

* $P_x$ is block-diagonal with $n_z$ identical $n_x \times n_x$ banded blocks (bandwidth $q_x$).

* $P_z$ is block-banded with bandwidth $q_z$. Each block is a scaled identity matrix.

* As a whole, $P_z$ and $P_{\lambda}$ may be viewed as banded matrices with bandwidth $q = q_z \times n_x$.

This implies that all statements made in the one-dimensional case carry over to the two-dimensional case with this value of $q$. It also suggests that, if $(q_x + 1)/n_x < (q_z + 1)/n_z$, dimensions $x$ and $z$ should be permuted before applying WH smoothing for maximal efficiency.

As in the one-dimensional case, the generalized Fellner-Schall update formula does not require the full computation of $V$. Indeed, to compute $\text{tr}[VP_x]$ and $\text{tr}[VP_z]$, we only need access to elements of $V$ for which either $P_x$ or $P_z$ is non-zero. From what precedes, $P_x$ has bandwidth $q_x$ while $P_z$ only contains $q_z$ non-zero diagonals on each side of the leading diagonal. As $V$ is symmetric, we only need to compute $q_x + q_z + 1$ diagonals of $V$ for an associated cost of $O([q_x + q_z]n^2)$ instead of $O(n^3)$.

With the Newton method, while computing $V = KK^T$ still incurs a $O(n^3)$ cost, matrix multiplications like $VP_j$ or $V(\partial W/\partial \rho_j + P_j)$ can be performed block-wise. It may easily be checked for example that the products $VP_x$ and $V(\partial W/\partial \rho_x + P_x)$ have a cost of $O(q_x n^2)$ while the products $VP_z$ and $V(\partial W/\partial \rho_z + P_z)$ have a cost of $O(q_z n^2)$.

## Summary of complexity gains

Thanks to the banded structure, most computations involved in WH smoothing can be accelerated by a factor of $n/(q+1)$ in the 1D case and $\max(n_x/(q_z+1), n_z/(q_x+1))$ in the 2D case. There are 3 notable exceptions:

* Cholesky decomposition is improved from $O(n^3)$ to $O(q^2n)$—a quadratic speed-up.

* Computation of $V = KK^T$ remains $O(n^3)$.

* Some matrix products required by Newton method get a full $n/(q_z+1)$ or $n/(q_x+1)$ speed-up in the 2D case.
The solution to the smoothing problem, as in the one-dimensional case, is given by:

$$ \hat{y} = U(U^T W U + S_\lambda)^{-1} U^T W y \quad \text{where} \quad S_\lambda = \lambda_x I_{n_z} \otimes \Sigma_x + \lambda_z \Sigma_z \otimes I_{n_x}. \quad (19) $$

Figure 8 represents the residual degrees of freedom associated with each parameter after applying the smoothing, in the two-dimensional case, for different combinations of the smoothing parameters. Similar to the one-dimensional case, these degrees of freedom decrease as the smoothing parameters increase and are particularly small for higher eigenvalues. The eigenvectors are sorted in ascending order of eigenvalues for each one-dimensional penalty matrix $D_{n_x,q_x}^T D_{n_x,q_x}$ and $D_{n_z,q_z}^T D_{n_z,q_z}$.

<figure><img src="image_10.png" /><figcaption>Figure 8: Distribution of residual degrees of freedom of the model by eigenvector after applying Whittaker-Henderson smoothing in the two-dimensional case, for different combinations of smoothing parameters.</figcaption></figure>
<table><caption>Table 3: Compared theoretical leading-order costs associated with the key steps in smoothing computations for several frameworks. All cells should be read as O(...).</caption><thead><tr><th>Computation</th><th>Dense</th><th>Banded</th><th>Rank-reduced</th></tr></thead><tbody><tr><td>X<sup>T</sup>WX</td><td>0</td><td>0</td><td>np<sup>2</sup></td></tr><tr><td>X<sup>T</sup>Wz</td><td>n</td><td>n</td><td>np</td></tr><tr><td>P<sub>λ</sub></td><td>n<sup>2</sup></td><td>qn</td><td>p</td></tr><tr><td>X<sup>T</sup>WX + P<sub>λ</sub></td><td>n</td><td>n</td><td>p</td></tr><tr><td>R</td><td>n<sup>3</sup></td><td>q<sup>2</sup>n</td><td>p<sup>3</sup></td></tr><tr><td>β̂<sub>λ</sub></td><td>n<sup>2</sup></td><td>qn</td><td>p<sup>2</sup></td></tr><tr><td>ŷ<sub>λ</sub> = Xβ̂<sub>λ</sub></td><td>0</td><td>0</td><td>np</td></tr><tr><td>ML/LAML</td><td>qn</td><td>qn</td><td>n</td></tr><tr><td>K = R<sup>-1</sup></td><td>n<sup>3</sup></td><td>qn<sup>2</sup></td><td>p<sup>3</sup></td></tr><tr><td>V = KK<sup>T</sup></td><td>n<sup>3</sup></td><td>n<sup>3</sup></td><td>p<sup>3</sup></td></tr><tr><td>Brent/Nelder-Mead</td><td>n<sup>3</sup></td><td>q<sup>2</sup>n</td><td>p<sup>3</sup></td></tr><tr><td>Fellner-Schall</td><td>n<sup>3</sup></td><td>qn<sup>2</sup></td><td>p<sup>3</sup></td></tr><tr><td>Newton</td><td>n<sup>3</sup></td><td>n<sup>3</sup></td><td>p<sup>3</sup></td></tr></tbody></table>

Table 3 summarises theoretical complexities across different frameworks, including a typical generalized additive model framework for which the penalization matrix is diagonal. This last framework is used by the rank-reduced WH smoothing approach introduced next, as well as the P-spline alternative used for comparison.

## Empirical gains

Figure 5 compares actual computation times of WH smoothing (two-dimensional, outer iteration), showing that adapting the implementation to exploit banded structures results in large speed gains:

* The Nelder-Mead method benefits the most, with a 25 × speedup compared to dense computation.

* Newton and Fellner-Schall methods see 6.6 × and 10 × improvements, respectively, making them fall behind the Nelder-Mead method.

As a final advantage, Brent and Nelder-Mead heuristic methods rely solely on banded matrices that can be stored as compact matrices of dimensions $(q + 1) \times n$, adding further efficiency.

## 6.4 Natural parameterization and rank reduction of WH smoothing

Demmler and Reinsch (1975) proposed a natural parameterization for penalized smoothers using the eigendecomposition of the penalization matrix. This provides both an intuitive
## 6.4 Eigendecomposition and rank reduction

In addition to providing a more intuitive interpretation of the smoothing effect, the eigendecomposition of the penalty matrix is key to reducing the dimension of the optimization problem associated with it. Figures 7 and 8 show that eigenvectors associated with higher eigenvalues are more penalized by the smoothing, to the extent that a large majority of the eigenvectors represented in these figures have residual degrees of freedom very close to 0 for intermediate values of $\lambda$. This suggests that if we simply remove the parameters associated with these eigenvectors from the model, the smoothing results would be minimally affected. We propose setting the coordinates associated with the $n-p$ largest eigenvectors to 0, thus retaining only a reduced number $p < n$ of parameters to estimate.

In the one-dimensional case, this translates to using instead of $\hat{y}$ the approximate estimator:

$$ \hat{y}_p = U_p (U_p^T W U_p + \lambda \Sigma_p)^{-1} U_p^T W y \quad (20) $$

where $U_p$ contains the first $p$ columns of $U$, and $\Sigma_p$ is a diagonal matrix constructed from the $p$ smallest eigenvalues of $\Sigma$. In the two-dimensional case, we use:

$$ \hat{y}_{p_x, p_z} = U_{p_x, p_z} (U_{p_x, p_z}^T W U_{p_x, p_z} + \lambda_x I_{p_z} \otimes \Sigma_{x, p_x} + \lambda_z \Sigma_{z, p_z} \otimes I_{p_z})^{-1} U_{p_x, p_z}^T W y \quad (21) $$

where $U_{p_x, p_z} = U_{z, p_z} \otimes U_{x, p_x}$, $U_{x, p_x}$ (resp. $U_{z, p_z}$) contains the first $p_x$ (resp. $p_z$) columns of $U_x$ (resp. $U_z$), and $\Sigma_{x, p_x}$ (resp. $\Sigma_{z, p_z}$) corresponds to a diagonal matrix constructed from the $p_x$ (resp. $p_z$) smallest eigenvalues of $\Sigma_x$ (resp. $\Sigma_z$).

In the two-dimensional case, there are several possible strategies for choosing the pair $(p_x, p_z)$. One simple solution is to set an upper bound $p_{max}$ on the number of parameters to be retained and choose $p_x$ and $p_z$ proportionally to the numbers $n_x$ and $n_z$ of categories in each dimension. Let $\kappa = \sqrt{p_{max}/n_x n_z}$, and choose $p_x = \lfloor \min(\kappa, 1)n_x \rfloor$ and $p_z = \lfloor \min(\kappa, 1)n_z \rfloor$. By construction, the pair $(p_x, p_z) \in \mathbb{N}^2$ defined in this way satisfies the condition $p_x \times p_z \le \kappa^2 n_x n_z = p_{max}$. Other strategies may be more effective. If it is known, for example, that the phenomenon under study exhibits greater smoothness in one dimension or the other, it may be interesting to retain a smaller proportion of components associated with that dimension.

The results presented here can be directly applied to generalized Whittaker-Henderson smoothing by replacing $y$ with $z_k$ and $W$ with $W_k$ in Equations 20 and 21. This approach allows for reducing the rank of the optimization problem to be solved and can be used at all stages of the algorithms presented. However, we limit its use to the selection of the smoothing parameter and prefer, once the optimal parameter is chosen, to use the estimator $\hat{y}$ (or $\hat{\theta}$) associated with the full-rank problem in Equation 3 or Algorithm 1.
<figure><img src="image_7.png" /><figcaption>Figure 5: Computation time comparison for 2D generalized WH smoothing with outer iteration. The speed-up factor is computed relative to the original dense method using the Nelder-Mead algorithm.</figcaption></figure>

interpretation of the smoothing mechanism and a foundation for dimension reduction via rank-restricted estimation.

## One-dimensional Case

In one dimension, let $D_{n,q}^T D_{n,q} = U \Sigma U^T$ be the eigendecomposition of the penalty matrix, where $U$ is orthogonal and $\Sigma$ diagonal with non-negative eigenvalues. A change of variable $\theta = U\beta$ transforms the WH optimization into:

$$ \hat{\mathbf{y}} = U\hat{\beta} \quad \text{where} \quad \hat{\beta} = \underset{\beta}{\operatorname{argmin}} \left\{ (\mathbf{y} - U\beta)^T W (\mathbf{y} - U\beta) + \lambda\beta^T \Sigma\beta \right\} $$

yielding the solution:

$$ \hat{\mathbf{y}} = U(U^T W U + S_{\lambda})^{-1} U^T W \mathbf{y} \quad \text{where} \quad S_{\lambda} = \lambda \Sigma. $$

This formulation shows that WH smoothing decomposes the signal into eigenvector components and attenuates each according to the associated eigenvalue—the higher the eigenvalue, the stronger the shrinkage.

We refer to Section E of the appendices for graphical illustrations of:

* the basis eigenvectors of $D_{n,q}^T D_{n,q}$;

* the evolution of their effective degrees of freedom under smoothing.

These figures show that only the first few components retain substantial degrees of freedom under moderate smoothing, motivating dimensionality reduction.
## 6.5 Impact of rank reduction method

Here, we aim to assess the impact of using the rank-reduced estimators defined in Section 6.4 on the smoothing results and computation time. To do so, we start from the 100 replicate datasets (one-dimensional and two-dimensional) used in Section 5.4. Since the use of rank reduction is limited to the selection of the smoothing parameter, and the final result is always obtained from the full-rank estimator, we can again use the error measure defined by Equation 14 to quantify the impact of the rank reduction method on the smoothing results. The quantity $\ell_{\text{ML}}^{m}(\lambda)$ involved in this formula will be computed based on the values of $\lambda$ obtained from the *performance iteration* approach for different values of the number of retained components $p$. The quantities $\ell_{\text{ML}}^{m}(\lambda_{\text{outer}})$ and $\ell_{\text{ML}}^{m}(\infty)$ appearing in the expression of the criterion will always be calculated using the *outer iteration* method and retaining all components in order to capture potential cross-effects of the choice of method and the application of rank reduction.

<figure><img src="image_11.png" /><figcaption>Figure 9: Distribution of the discrepancies induced by the proposed rank reduction method, depending on the number of retained eigenvectors during the selection of the smoothing parameter, in the one-dimensional and two-dimensional case.</figcaption></figure>

Figure 9 represents the impact of rank reduction on the smoothing results as a function of the number of retained components. In the one-dimensional case, reducing the number of retained eigenvectors from 74 to 10 has only a marginal impact on the results obtained. However, as shown in Figure 10, it reduces the computation time by a factor of 3. In the two-dimensional case, reducing the number of components to 160 reduces the computation time by a factor close to 10, with imperceptible differences in the obtained results. The proposed rank reduction method, in combination with the *performance iteration* approach, thus allows, in the two-dimensional case, by considering only 160 components, to reduce the computation time by a factor of 40 at the cost of an average error of 0.03%, which exceeds 0.5% in only one of the 100 simulated datasets considered. This solution offers significant operational advantages, especially for applications mentioned in Section 6.1 where computation times can be much
## Two-dimensional case

In two dimensions, the penalization matrix takes the form:

$$P_{\lambda} = \lambda_{x}I_{n_{x}} \otimes D_{n_{x},q_{x}}^{T}D_{n_{x},q_{x}} + \lambda_{z}D_{n_{z},q_{z}}^{T}D_{n_{z},q_{z}} \otimes I_{n_{x}},$$

with eigendecompositions $D_{n_x,q_x}^T D_{n_x,q_x} = U_x \Sigma_x U_x^T$ and $D_{n_z,q_z}^T D_{n_z,q_z} = U_z \Sigma_z U_z^T$.

Define $U = U_z \otimes U_x$, and $\theta = U\beta$. Then the WH estimate becomes:

$$\hat{\mathbf{y}} = U(U^T W U + S_{\lambda})^{-1} U^T W \mathbf{y} \quad \text{where} \quad S_{\lambda} = \lambda_x I_{n_z} \otimes \Sigma_x + \lambda_z \Sigma_z \otimes I_{n_x}.$$

As in the one-dimensional case, this representation reveals how smoothing operates via coordinate-wise shrinkage in the eigenbasis. Section E of the appendices displays the corresponding per-parameter effective degrees of freedom.

## Rank reduction strategy

Inspection of the effective degrees of freedom reveals that many components are heavily shrunk, especially those associated with high eigenvalues. This suggests reducing the dimension by keeping only the $p < n$ components with the lowest eigenvalues.

In the one-dimensional case, the reduced-rank approximation is:

$$\hat{\mathbf{y}}_p = U_p(U_p^T W U_p + \lambda \Sigma_p)^{-1} U_p^T W \mathbf{y}$$

where $U_p$ and $\Sigma_p$ consist of the first $p$ eigenvectors and their corresponding eigenvalues respectively.

In the two-dimensional case, we retain $p_x$ and $p_z$ eigenvectors in each dimension and use:

$$\hat{\mathbf{y}}_{p_x, p_z} = U_{p_x, p_z}(U_{p_x, p_z}^T W U_{p_x, p_z} + \lambda_x I_{p_z} \otimes \Sigma_{x, p_x} + \lambda_z \Sigma_{z, p_z} \otimes I_{p_x})^{-1} U_{p_x, p_z}^T W \mathbf{y}$$

with $U_{p_x, p_z} = U_{z, p_z} \otimes U_{x, p_x}$. In that case, given a target number of parameters $p_{\max}$, we propose selecting $(p_x, p_z)$ such that $p_x p_z \le p_{\max}$ and $p_x/n_x \approx p_z/n_z$ using the rule:

$$\kappa = \sqrt{p_{\max}/n_x n_z}, \quad p_x = \lfloor \min(\kappa, 1)n_x \rfloor, \quad p_z = \lfloor \min(\kappa, 1)n_z \rfloor.$$

Adaptations for generalized WH smoothing follow by replacing $(\mathbf{y}, W)$ with $(\mathbf{z}_k, W_k)$ in the above expressions.
more significant than in the analysis presented here.

<figure><img src="image_12.png" /><figcaption>Figure 10: Distribution of computation time for the proposed rank reduction method, depending on the number of retained eigenvectors during the selection of the smoothing parameter, in the one-dimensional and two-dimensional case.</figcaption></figure>

# 7 How to extrapolate the smoothing?

## 7.1 Defining the (unconstrained) extrapolation of the smoothing

Let $\hat{y}$ be the result of Whittaker-Henderson smoothing for an observation vector $y$, which corresponds in the one-dimensional case to an explanatory variable vector $x$ or in the two-dimensional case to a combination of explanatory variables $x$ and $z$. Suppose we have an observation vector $x_+$ or combinations of vectors $x_+$ and $z_+$ for which we want to make predictions using the model. Since Whittaker-Henderson smoothing applies only to evenly spaced observations, without loss of generality, let $x_+$ and $z_+$ be sequences of consecutive integers such that $x \subset x_+$ and $z \subset z_+$. In the one-dimensional case, let $n_+$ be the length of $x_+$, and in the two-dimensional case, let $n_{x+}$ and $n_{z+}$ be the respective lengths of $x_+$ and $z_+$, with $n_+ = n_{x+} \times n_{z+}$. Let $C_x$ (resp. $C_z$) be a matrix of size $n_x \times n_{x+}$ (resp. $n_z \times n_{z+}$) defined as $C_j = (0_{n_{j,\min j-\min j_+}|I_{n_{j,\min j-\max j_+}}|0_{n_{j,\max j_+-\max j}}}$ for $j \in \{x, z\}$, and let us note:

$$C = \begin{cases} C_x & \text{in the one-dimensional case} \\ C_z \otimes C_x & \text{in the two-dimensional case} \end{cases}$$
## Efficient computation via GLAM

Currie, Durban, and Eilers (2006) propose a general framework, Generalized Linear Array Models (GLAM), that exploits Kronecker structure for efficient computations. In our context, the model matrix $U_{p_x, p_z}$ inherits a Kronecker product form, allowing operations that rely on this matrix to be executed dimension-wise without explicit construction of the full matrix. This significantly reduces memory use and computation time in the two-dimensional rank-reduced WH framework.

## Impact of using the rank-reduced basis

We now evaluate the impact of the rank-reduced WH basis introduced in Section 6.4 in terms of both smoothing accuracy and computational speed. For context, results are compared against those obtained using P-spline smoothing with the same number of basis functions.

To ensure a fair comparison, both approaches were implemented in the same computational framework, including the use of GLAM in the two-dimensional case—only the structure of the basis (and hence the model matrix) differs. The penalty structure, as well as the unpenalized fixed effects (polynomials of degree $q-1$), are identical.

In addition to the full basis of size 450 ($30 \times 15$), three reduced basis of respective size 288 ($24 \times 12$), 128 ($16 \times 8$) and 32 ($8 \times 4$) were considered.

As in Section 5.6, accuracy is assessed using the relative LAML error defined in Equation 13. Note, however, that since both reduced-rank and P-spline smoothers rely on different bases and penalization matrices, their LAML expressions are different from the one used for full-rank WH smoothing. Hence, a reduced model can exhibit a higher LAML than the full-rank version at its selected smoothing parameter.

Figure 6 summarises the average speed-up achieved by both the reduced-rank and P-spline smoothers compared to the full-rank WH smoothing. As the number of retained parameters decreases, computation time drops substantially. Compared to the full-rank WH smoothing (unoptimized):

* the 128-parameter basis achieves an 88 × speed-up;

* the 32-parameter basis achieves up to 256 × faster computation.

The alternate iteration and performance iteration strategies outperform the outer iteration in the reduced setting, primarily because model matrix construction becomes the new computational bottleneck—even with the use of the GLAM framework. In this context, the Newton algorithm combined with alternate iteration proves to be the most efficient, with the generalized Fellner-Schall update being nearly as competitive for smaller bases.

The gains in computational speed come with a moderate tradeoff in estimation accuracy. As shown in Figure 7, the relative LAML errors remain small:
The matrix C of size n × n+ defined in this way has the following properties:

* For any vector $\mathbf{y}_+$ of size $n_+$, $C\mathbf{y}_+$ is a vector of size n containing only the n values of $\mathbf{y}_+$ matching the positions of the observations in the initial smoothing,

* $C^T\mathbf{y}$ is a vector of size $n_+$ containing the values of $\mathbf{y}$ at the positions corresponding to the initial observations and zeros elsewhere,

* $CC^T = I_n$.

The extrapolation of Whittaker-Henderson smoothing can now be defined as finding the solution $\hat{\mathbf{y}}_+$ to the extended optimization problem:

$$ \hat{\mathbf{y}}_+ = \underset{\boldsymbol{\theta}_+}{\operatorname{argmin}} \left\{ (\mathbf{y}_+ - \boldsymbol{\theta}_+)^T W_+ (\mathbf{y}_+ - \boldsymbol{\theta}_+) + \boldsymbol{\theta}_+^T P_+ \boldsymbol{\theta}_+ \right\} \quad (22) $$

where $\mathbf{y}_+ = C^T\mathbf{y}$, $W_+ = C^TWC$, and

$$ P_+ = \begin{cases} \lambda D_{n+,q}^T D_{n+,q} & \text{in the one-dimensional case} \\ \lambda_x I_{z+} \otimes D_{n_{x+},q_x}^T D_{n_{x+},q_x} + \lambda_z D_{n_{z+},q_z}^T D_{n_{z+},q_z} \otimes I_{x+} & \text{in the two-dimensional case} \end{cases} $$

In this last expression, $\lambda$, $\lambda_x$ and $\lambda_z$ are the parameters chosen for the smoothing. Similar to the initial smoothing, Equation 22 involves a fidelity criterion and a smoothness criterion. It should be noted that:

$$ (\mathbf{y}_+ - \boldsymbol{\theta}_+)^T W_+ (\mathbf{y}_+ - \boldsymbol{\theta}_+) = (C^T\mathbf{y} - \boldsymbol{\theta}_+)^T C^TWC (C^T\mathbf{y} - \boldsymbol{\theta}_+) = (\mathbf{y} - \boldsymbol{\theta})^T W (\mathbf{y} - \boldsymbol{\theta}) $$

hence the fidelity criterion remains unchanged compared to the initial estimation problem. This is consistent with the fact that extrapolation does not involve any additional data. The smoothness criterion, however, applies to all elements of the vector $\boldsymbol{\theta}_+$, whether initially present or not.

## 7.2 Unconstrained solution for the one-dimensional case

The solution to the extended optimization problem of Equation 22 is readily obtained by taking the derivatives in $\boldsymbol{\theta}_+$ and setting them to 0, as in Section 1.5. This yields the solution:

$$ \hat{\mathbf{y}}_+ = (W_+ + P_+)^{-1} W_+ \mathbf{y}_+ \quad \text{where} \quad \mathbf{y}_+ = C^T\mathbf{y} \quad \text{and} \quad W_+ = C^TWC \quad (23) $$

Let us further assume that $\mathbf{y}_+|\boldsymbol{\theta}_+ \sim N(\boldsymbol{\theta}_+, W_+)$ where $\boldsymbol{\theta}_+ \sim N(0, P_+)$ and proceed as in Section 2 to obtain the following credible intervals:
<figure><img src="image_8.png" /><figcaption>Figure 6: Computation speed improvement from WH smoothing with a reduced-rank basis (solid lines) or P-spline basis (dashed lines), relative to unoptimized full-rank WH smoothing, as a function of basis size.</figcaption></figure>

* For the 128-parameter basis, the average error is just 0.82%.

* For the 32-parameter basis, it rises to 2.26%.

Across all sizes, the reduced-rank WH smoother slightly outperforms the P-spline smoother in terms of LAML error, confirming its effectiveness as a principled dimension reduction strategy.
$$ \mathbb{E}(\mathbf{y}_+) | \mathbf{y}_+ \in \left[ (W_+ + P_+)^{-1} W_+ \mathbf{y}_+ \pm \Phi^{-1} \left(1 - \frac{\alpha}{2}\right) \sqrt{\mathrm{diag} \{(W_+ + P_+)^{-1}\}} \right] \quad (24) $$

To get a better understanding about how the variance-covariance matrix $\Psi_+ = (W_+ + P_+)^{-1}$ for the unconstrained extrapolation problem of Equation 22 is related to the variance-covariance matrix $\Psi = (W + P_\lambda)^{-1}$ of the original smoothing problem, let us introduce:

$$ \bar{C}_j = \begin{bmatrix} I_{\min j - \min j_+, \min j - \min j_+} & 0_{\min j - \min j_+, n_j} & 0_{\min j - \min j_+, \max j_+ - \max j} \\ 0_{\max j_+ - \max j, \min j - \min j_+} & 0_{\max j_+ - \max j, n_j} & I_{\max j_+ - \max j, \max j_+ - \max j} \end{bmatrix} \quad \text{for } j \in \{x, z\}, $$

$$ \bar{C} = \begin{cases} \bar{C}_x & \text{in the one-dimensional case} \\ \bar{C}_z \otimes \bar{C}_x & \text{in the two-dimensional case} \end{cases} \quad \text{and} \quad Q = \begin{bmatrix} C \\ \bar{C} \end{bmatrix}. $$

With this definition, $Q$ is a permutation matrix of size $n_+ \times n_+$ which select the rows whose indices correspond to those of the initial observation positions in the extended observation vector and put them first. It also verifies $Q^{-1} = Q^T = Q$.

In the unidimensional case, the extended difference matrix $D_{n_+,q}$ takes the form:

$$ D_{n_+,q} = \begin{bmatrix} D_{2-} & D_{1-} & 0 \\ 0 & D_{n,q} & 0 \\ 0 & D_{1+} & D_{2+} \end{bmatrix} = Q \begin{bmatrix} D_{n,q} & 0 & 0 \\ D_{1-} & D_{2-} & 0 \\ D_{1+} & 0 & D_{2+} \end{bmatrix} Q = Q \begin{bmatrix} D_{n,q} & 0 \\ D_{1} & D_{2} \end{bmatrix} Q $$

$$ \text{where } D_1 = \begin{bmatrix} D_{1-} \\ D_{1+} \end{bmatrix} \quad \text{and} \quad D_2 = \begin{bmatrix} D_{2-} & 0 \\ 0 & D_{2+} \end{bmatrix} $$

Detailed expressions of matrices $D_{1+}$ and $D_{2+}$ for the most common values of $q$ may be found in Carballo et al. (2021) where they are simply noted $D_1$ and $D_2$. To cover extrapolation on both sides of the initial observation vector, matrices $D_{1-}$ and $D_{2-}$ were introduced. Those may be simply obtained by taking the transpose of adequate size $D_1$ and $D_2$ matrices found in Carballo et al. (2021). The extended weight and penalization matrices thus take the form:

$$ W_+ = Q \begin{bmatrix} W & 0 \\ 0 & 0 \end{bmatrix} Q \quad \text{and} \quad P_+ = D_{n_+,q}^T D_{n_+,q} = \lambda Q \begin{bmatrix} P_\lambda + P_+^{11} & P_+^{12} \\ P_+^{21} & P_+^{22} \end{bmatrix} Q \quad (25) $$

$$ \text{where } P_+^{11} = \lambda D_1^T D_1, \quad P_+^{12} = \lambda D_1^T D_2, \quad P_+^{21} = \lambda D_2^T D_1 \text{ and } P_+^{22} = \lambda D_2^T D_2. $$

Applying the formulas for the inversion of a symmetric matrix partitioned with 2 × 2 blocks to $Q(W_+ + P_+)Q$, we obtain a more detailed expression for $\Psi_+$:
<figure><img src="image_9.png" /><figcaption>Figure 7: Relative LAML error of WH smoothing with a reduced-rank basis (solid lines) or a P-spline basis (dashed lines), with respect to unoptimized full-rank WH smoothing, as a function of basis size.</figcaption></figure>

# 7 How to extrapolate the smoothing?

Semi-parametric methods such as P-splines and Whittaker-Henderson (WH) smoothing naturally allow for extrapolation—that is, predicting values outside the range of the original data. Extrapolation is handled by solving an extended smoothing problem where extrapolated positions are associated with zero-weight observations.

However, in the two-dimensional case, extrapolation must be performed carefully: constraints are needed to ensure that the extrapolated solution remains consistent with the original smoothing result over the observed data. Following the approach introduced by Carballo, Durban, and Lee (2021) for P-splines, we now extend WH smoothing to support extrapolation while also enabling the construction of credibility intervals that capture uncertainty both inside and outside the original observation domain.

## 7.1 Defining the extrapolation of the smoothing

Let $\hat{y}$ be the WH smoothing result obtained from an observation vector $y$ defined over positions $x$ (in 1D) or $(x, z)$ (in 2D). We wish to extend predictions to a larger domain $x^+$ (or $(x^+, z^+)$), with $x \subset x^+$ and similarly for $z$.

To preserve WH smoothing's requirement for evenly spaced points, we assume that $x^+$ and $z^+$ are sequences of consecutive integers. Let $n_+$ be the length of $x^+$ in the on-dimensional
$$ \Psi_+ = Q \begin{bmatrix} \Psi_+^{11} & -\Psi_+^{11} P_+^{12} (P_+^{22})^{-1} \\ -(P_+^{22})^{-1} P_+^{21} \Psi_+^{11} & (P_+^{22})^{-1} + (P_+^{22})^{-1} P_+^{21} \Psi_+^{11} P_+^{12} (P_+^{22})^{-1} \end{bmatrix} Q. \quad (26) $$

where $\Psi_+^{11} = [W + P_\lambda + P_+^{11} - P_+^{12}(P_+^{22})^{-1}P_+^{21}]^{-1}$.

Let us denote by initial positions coefficients the subvector $C\hat{\mathbf{y}}_+$ that correspond to the coefficients at the position of the initial observations in $\hat{\mathbf{y}}_+$ and by new positions coefficients the remaining coefficients $\bar{C}\hat{\mathbf{y}}_+$. The initial positions coefficients may be recovered as $C\hat{\mathbf{y}}_+ = C\Psi_+ C^T W_+ y_+ = \Psi_+^{11} Wy$. This does not simplify to $\Psi Wy = \hat{\mathbf{y}}$ for all weighted observation vectors $Wy$ unless $\Psi_+^{11} = \Psi$. Indeed, the initial positions coefficients are chosen to optimize the overall smoothness of the extrapolated coefficient vectors and note just the smoothness of $C\hat{\mathbf{y}}_+$. Besides, the expression of $\Psi_+^{22}$, it contains two terms: an innovation error $(P_+^{22})^{-1}$ associated with the prior on the new positions coefficients while $(P_+^{22})^{-1}P_+^{21}\Psi_+^{11}P_+^{12}(P_+^{22})^{-1}$ represents an additional uncertainty on the new positions coefficients caused by the uncertainty on the initial positions coefficients.

In the one-dimensional case, $D_2$ is a block-diagonal matrix of triangular matrices with non-zero diagonal elements and is thus non-singular. Hence $P_+^{11} - P_+^{12}(P_+^{22})^{-1}P_+^{21} = D_1^T D_1 - D_1^T D_2(D_2^T D_2)^{-1}D_2^T D_1 = 0$ which means that $\Psi_+^{11} = (W + P_\lambda)^{-1} = \Psi$ and $C\hat{\mathbf{y}}_+ = \hat{\mathbf{y}}$. Therefore, in the one-dimensional case, the solution given by Equation 22 preserves the values from the original fit. As shown by Carballo et al. (2021), it is indeed always possible in the case of penalizations based on difference matrices to pick the new positions coefficients so that the smoothness criterion does not increase or in other words to find a perfectly smooth extrapolation for the fit. Figure 11 shows the extrapolation associated with the Whittaker-Henderson smoothing applied to the data use to produce Figure 1.

## 7.3 Constrained solution for the two-dimensional case

In the two-dimensional case, while the extended penalization matrix $P_+$ still takes the form of Equation 25, expressions of $P_+^{11}$, $P_+^{12}$, $P_+^{21}$ and $P_+^{22}$ are more complex. In particular, $P_+^{11} - P_+^{12}(P_+^{22})^{-1}P_+^{21} \neq 0$ which implies that $\Psi_+^{11} \neq \Psi$ and $C\hat{\mathbf{y}}_+ \neq \hat{\mathbf{y}}$. Solving Equation 23 thus leads to a change in the value of the initial positions coefficients compared to the coefficients obtained during the initial smoothing as shown by Carballo, Durban, and Lee (2021). Indeed, the smoothness criterion includes penalizations on both rows and columns and it is no longer possible, as in the one-dimensional case, to extrapolate the fit without increasing that criterion. As the smoothness criterion carries more weight in the extrapolation problem compared to the original problem, the optimal solution to the extended optimization problem will compromise on the fidelity to the initial observations in order to improve the overall smoothness.

To obtain an estimator $\hat{\mathbf{y}}_+^*$ that minimizes the penalized regression problem under the constraint of preserving the initial coefficients, i.e. $C\hat{\mathbf{y}}_+^* = \hat{\mathbf{y}}$, we follow the approach proposed by Carballo,
case. In the two-dimensional case, let $n_{x+}$ and $n_{z+}$ be the lengths of $\mathbf{x}+$ and $\mathbf{z}+$ and note $n_+ = n_{x+} \times n_{z+}$.

We define matrices $C_x$ and $C_z$ such that each extracts the indices of the original data from the larger domain. Specifically: $C_j = (O | I_{n_j} | O)$, where $j \in \{x, z\}$ and $I_{n_j}$ is an identity matrix aligned with the observed positions. Define the matrix $C$ as:

$$C = \begin{cases} C_x & \text{in the one-dimensional case,} \\ C_z \otimes C_x & \text{in the two-dimensional case.} \end{cases}$$

Then $C$ has the following useful properties:

* For any full-domain vector $\mathbf{y}+$, $C\mathbf{y}+$ returns the observed values only.

* $C^T\mathbf{y}$ embeds the observed values into a larger zero-padded vector.

* $CC^T = I_n$ and $C^T C$ is a $2 \times 2$ block matrix with an identity matrix block and zeros everywhere else.

The extrapolated WH smoothing is defined as the solution to the following extended problem:

$$\hat{\mathbf{y}}_+ = \underset{\mathbf{\theta}_+}{\operatorname{argmin}} \left\{ (\mathbf{y}_+ - \mathbf{\theta}_+)^T W_+ (\mathbf{y}_+ - \mathbf{\theta}_+) + \mathbf{\theta}_+^T P_+ \mathbf{\theta}_+ \right\} \quad (14)$$

where:

* $\mathbf{y}_+ = C^T\mathbf{y}$ is the extended data vector (zeros for unobserved points),

* $W_+ = C^T W C$ is the extended weight matrix (zeros for unobserved points),

* $P_+$ is the penalization matrix over the extended grid, defined as:

$$P_+ = \begin{cases} \lambda D_{n_+,q}^T D_{n_+,q} & \text{in the one-dimensional case,} \\ \lambda_x I_{z+} \otimes D_{n_{x+},q_x}^T D_{n_{x+},q_x} + \lambda_z D_{n_{z+},q_z}^T D_{n_{z+},q_z} \otimes I_{x+} & \text{in the two-dimensional case.} \end{cases}$$

Importantly, the smoothing parameters $\lambda$, $\lambda_x$, and $\lambda_z$ must remain fixed during extrapolation— they are inherited from the original fit and no new information is introduced.

The fidelity term in Equation 14 simplifies to:

$$(\mathbf{y}_+ - \mathbf{\theta}_+)^T W_+ (\mathbf{y}_+ - \mathbf{\theta}_+) = (C^T \mathbf{y} - \mathbf{\theta}_+)^T C^T W C (C^T \mathbf{y} - \mathbf{\theta}_+) = (\mathbf{y} - \mathbf{\theta})^T W (\mathbf{y} - \mathbf{\theta})$$

where $\mathbf{\theta} = C\mathbf{\theta}_+$. This is the fidelity term from the original fit.

The smoothness criterion, on the other hand, now applies to the entire extended domain, constraining the extrapolated parts of $\hat{\mathbf{y}}_+$ to remain smooth and consistent with the trend learned from the data.
<figure><img src="image_13.png" /><figcaption>Figure 11: Extrapolation of Whittaker-Henderson smoothing in the one-dimensional case. The smoother extrapolation is achieved, on both sides of the smoothing, using a polynomial of degree q - 1 (in this case a straight line as q = 2).</figcaption></figure>

Durban, and Lee (2021) and introduce the Lagrange multiplier $\omega$. The associated constrained extended optimization problem is now written as:

$$ (\hat{\mathbf{y}}_+^*, \hat{\omega}) = \underset{\mathbf{\theta}_+^*, \omega}{\operatorname{argmin}} \left\{ (\mathbf{y}_+ - \mathbf{\theta}_+^*)^T W_+(\mathbf{y}_+ - \mathbf{\theta}_+^*) + \mathbf{\theta}_+^{*T} P_+ \mathbf{\theta}_+^* + 2\omega^T (C \mathbf{\theta}_+^* - \hat{\mathbf{y}}) \right\}. \quad (27) $$

Taking the partial derivatives of Equation 27 with respect to $\mathbf{\theta}_+^*$ and $\omega$ gives:

$$ \frac{\partial}{\partial \mathbf{\theta}_+^*} \left\{ (\mathbf{y}_+ - \mathbf{\theta}_+^*)^T W_+(\mathbf{y}_+ - \mathbf{\theta}_+^*) + \mathbf{\theta}_+^{*T} P_+ \mathbf{\theta}_+^* + 2\omega^T (C \mathbf{\theta}_+^* - \hat{\mathbf{y}}) \right\} = -2W_+(\mathbf{y}_+ - \mathbf{\theta}_+^*) + 2P_+ \mathbf{\theta}_+^* + 2\omega^T C $$

$$ \frac{\partial}{\partial \omega} \left\{ (\mathbf{y}_+ - \mathbf{\theta}_+^*)^T W_+(\mathbf{y}_+ - \mathbf{\theta}_+^*) + \mathbf{\theta}_+^{*T} P_+ \mathbf{\theta}_+^* + 2\omega^T (C \mathbf{\theta}_+^* - \hat{\mathbf{y}}) \right\} = 2(C \mathbf{\theta}_+^* - \hat{\mathbf{y}}) $$

Setting these derivatives to zero yields the linear system:

$$ \begin{bmatrix} W_+ + P_+ & C^T \\ C & 0 \end{bmatrix} \begin{bmatrix} \hat{\mathbf{y}}_+^* \\ \hat{\omega} \end{bmatrix} = \begin{bmatrix} W_+ \mathbf{y}_+ \\ \hat{\mathbf{y}} \end{bmatrix} $$

The solution for $\hat{\mathbf{y}}_+^*$ can be derived using formulas for the inversion of a symmetric partitioned matrix with 2 × 2 blocks:
The same extrapolation approach applies directly to generalized WH smoothing, simply by replacing $\mathbf{y}$ by $\mathbf{z}_k$ and $W$ by $W_k$, obtained at convergence of the PIRLS algorithm and setting $\sigma^2 = 1$ in the derived credible intervals.

## 7.2 Unconstrained solution for the 1D case

The solution to the extrapolation problem in Equation 14 can be obtained directly, as in Section 1.1, by taking derivatives with respect to $\theta_+$ and setting them to zero. This yields the closed-form solution:

$$ \hat{\mathbf{y}}_+ = (W_+ + P_+)^{-1} W_+ \mathbf{y}_+ \quad \text{where} \quad \mathbf{y}_+ = C^T \mathbf{y} \quad \text{and} \quad W_+ = C^T W C. $$

Assuming a Bayesian model where $\mathbf{y}+ | \theta_+ \sim \mathcal{N}(\theta_+, \sigma^2 W + \mathbf{1}^{-1})$ and $\theta_+ \sim \mathcal{N}(0, \sigma^2 P + \mathbf{1}^{-1})$, we obtain, as in Section 2, the following credible interval:

$$ \mathbb{E}(\mathbf{y}_+) | \mathbf{y}_+ \in \left[ (W_+ + P_+)^{-1} W_+ \mathbf{y}_+ \pm \Phi^{-1} \left(1 - \alpha/2\right) \sqrt{\sigma^2 \mathbf{diag}\{(W_+ + P_+)^{-1}\}} \right]. $$

To get a better understanding about how the variance-covariance matrix $V_+ = (W_+ + P_+)^{-1}$ for the unconstrained extrapolation problem of Equation 14 is related to the variance-covariance matrix $V = (W + P_\lambda)^{-1}$ of the original smoothing problem, introduce matrices $\bar{C}_j$ (for $j \in \{x, z\}$) which selects the rows in the extrapolated domain that are not part of the original data and define:

$$ \bar{C} = \begin{cases} \bar{C}_x & \text{in the one-dimensional case,} \\ \bar{C}_z \otimes \bar{C}_x & \text{in the two-dimensional case,} \end{cases} \quad \text{and} \quad Q = \left[ \frac{C}{\bar{C}} \right]. $$

With this definition, $Q$ is a permutation matrix moving observed positions to the top.

In the unidimensional case, the extended difference matrix $D_{n+,q}$ takes the block-wise form:

$$ D_{n+,q} = \begin{bmatrix} D_{2-} & D_{1-} & 0 \\ 0 & D_{n,q} & 0 \\ 0 & D_{1+} & D_{2+} \end{bmatrix} = Q^T \begin{bmatrix} D_{n,q} & 0 \\ D_1 & D_2 \end{bmatrix} Q \quad \text{where} \quad D_1 = \begin{bmatrix} D_{1-} \\ D_{1+} \end{bmatrix} \quad \text{and} \quad D_2 = \begin{bmatrix} D_{2-} & 0 \\ 0 & D_{2+} \end{bmatrix}. $$

The extended weight and penalization matrices may be rewritten:

$$ W_+ = Q^T \begin{bmatrix} W & 0 \\ 0 & 0 \end{bmatrix} Q \quad \text{and} \quad P_+ = D_{n+,q}^T D_{n+,q} = \lambda Q^T \begin{bmatrix} P_\lambda + P_+^{11} & P_+^{12} \\ P_+^{21} & P_+^{22} \end{bmatrix} Q $$

where $P_+^{ij} = \lambda D_i^T D_j$, for $i, j \in \{1, 2\}$.

This block structure allows us to apply standard results for partitioned matrix inverses to derive:
$$ \hat{\mathbf{y}}_{+}^{*} = (W_{+} + P_{+})^{-1} \left\{ I - C^{T}[C(W_{+} + P_{+})^{-1}C^{T}]^{-1}C(W_{+} + P_{+})^{-1} \right\} W_{+} \mathbf{y}_{+} + (W_{+} + P_{+})^{-1} C^{T}[C(W_{+} + P_{+})^{-1}C^{T}]^{-1} \hat{\mathbf{y}} $$

Since $W_+ = C^T W C$, the first term is actually zero, and this expression simplifies to:

$$ \hat{\mathbf{y}}_{+}^{*} = (W_{+} + P_{+})^{-1} C^{T} [C(W_{+} + P_{+})^{-1} C^{T}]^{-1} \hat{\mathbf{y}} = Q \left[ \Psi^{21} (\Psi^{11})^{-1} \right] \hat{\mathbf{y}} = Q \left[ -(P_{+}^{22})^{-1} P_{+}^{21} \right] \hat{\mathbf{y}} \quad (28) $$

which is a linear transformation of $\hat{\mathbf{y}}$. Defining $A_+^* = Q \begin{bmatrix} I \\ -(P_+^{22})^{-1} P_+^{21} \end{bmatrix}$, a natural candidate for the variance-covariance of $\mathbf{y}_+^*|\theta_+$ is given by:

$$ A_+^* \Psi A_+^{*T} = Q \begin{bmatrix} \Psi & -\Psi P_+^{12} (P_+^{22})^{-1} \\ -(P_+^{22})^{-1} P_+^{21} \Psi & (P_+^{22})^{-1} P_+^{21} \Psi P_+^{12} (P_+^{22})^{-1} \end{bmatrix} Q \quad (29) $$

Equation 29 is very similar to Equation 26 with however two differences. First, every occurrence of $\Psi_+^{11}$ is replaced by $\Psi$. This is consistent with the constraint that the initial positions coefficients are forced to take the value of the coefficients used during the initial smoothing. Second, as the solution to the constrained extended optimization problem of Equation 27 was expressed as a linear transformation of $\hat{\mathbf{y}}$, Equation 29 is missing the innovation error term $(P_+^{22})^{-1}$ associated with the prior on the new positions coefficients. Not including this term would be tantamount to considering that $\theta_+$ has some degree of wigglyness in the region of the initial data but is perfectly smooth everywhere else. Adding the innovation error back, we obtain the following variance-covariance matrix for the constrained optimization problem:

$$ \Psi_+^* = Q \begin{bmatrix} \Psi & -\Psi P_+^{12} (P_+^{22})^{-1} \\ -(P_+^{22})^{-1} P_+^{21} \Psi & (P_+^{22})^{-1} + (P_+^{22})^{-1} P_+^{21} \Psi P_+^{12} (P_+^{22})^{-1} \end{bmatrix} Q \quad (30) $$

which still verifies

$$ \Psi_+^* W_+ \mathbf{y}_+ = \Psi_+^* C W \mathbf{y} = Q \begin{bmatrix} I \\ -(P_+^{22})^{-1} P_+^{21} \end{bmatrix} \Psi W \mathbf{y} = \hat{\mathbf{y}}_+^* $$

The associated credible intervals are readily obtained as:

$$ \mathbb{E}(\mathbf{y}_+) | \mathbf{y}_+ \in \left[ \hat{\mathbf{y}}_+^* \pm \Phi^{-1} \left( 1 - \frac{\alpha}{2} \right) \sqrt{\mathrm{diag}(\Psi_+^*)} \right] $$
$$V_+ = Q^T \begin{bmatrix} V_+^{11} & V_+^{12} \\ V_+^{21} & V_+^{22} \end{bmatrix} Q = Q^T \begin{bmatrix} V_+^{11} & -V_+^{11} P_+^{12} (P_+^{22})^{-1} \\ -(P_+^{22})^{-1} P_+^{21} V_+^{11} & (P_+^{22})^{-1} P_+^{21} V_+^{11} P_+^{12} (P_+^{22})^{-1} + (P_+^{22})^{-1} \end{bmatrix} Q$$

with $V_+^{11} = [W + P_\lambda + P_+^{11} - P_+^{12}(P_+^{22})^{-1}P_+^{21}]^{-1}$.

From the above, we retrieve:

$$C\hat{y}_+ = CV_+W_+y_+ = CQ^T V_+QC^T Wy = V_+^{11}Wy.$$

This coincides with the original fit $\hat{y}$ only if $V_+^{11} = V$. In general, this equality does not hold, since the extrapolation solution minimizes the total smoothness of the extended vector, not just of the observed part.

In $V_+^{22}$, we identify:

* a propagation term: $(P_+^{22})^{-1}P_+^{21}V_+^{11}P_+^{12}(P_+^{22})^{-1}$, capturing the uncertainty transferred from the known part to the extrapolated part;

* an innovation error term: $(P_+^{22})^{-1}$ associated with the prior on the extrapolated coefficients $\bar{C}\hat{y}_+$.

In the one-dimensional case, $D_2$ is block-diagonal with invertible triangular blocks, so:

$$P_+^{11} - P_+^{12}(P_+^{22})^{-1}P_+^{21} = D_1^T D_1 - D_1^T D_2(D_2^T D_2)^{-1}D_2^T D_1 = 0$$

which means that $V_+^{11} = (W + P_\lambda)^{-1} = V$. This confirms the result from Carballo et al. (2021), namely that with a difference-based penalty, a perfectly smooth extrapolation that leaves the original fit unchanged can always be constructed in the one-dimensional case.

This behaviour is illustrated in Figure 8, which shows the extrapolated fit (with $q = 2$) obtained from generalized WH smoothing applied to the annuity portfolio used previously. The extrapolation follows a straight line—the polynomial of degree $q - 1 = 1$—and joins smoothly with the original curve.

## 7.3 Constrained solution for the 2D case

In the two-dimensional case, while the extended penalization matrix $P_+$ still takes the same structure as previously described, the expressions of its block components $P_+^{11}$, $P_+^{12}$, $P_+^{21}$ and $P_+^{22}$ are more complex. In particular, we no longer have the simplification $P_+^{11} - P_+^{12}(P_+^{22})^{-1}P_+^{21} = 0$, therefore $V_+^{11} \neq V$ and $C\hat{y}_+ \neq \hat{y}$. Solving the unconstrained extrapolation problem thus leads to a modification of the estimated coefficients for the observed data positions, as demonstrated by Carballo, Durban, and Lee (2021).

This difference arises because, unlike the one-dimensional case, the smoothness criterion in two dimensions penalizes both rows and columns simultaneously, making it impossible to extrapolate without increasing the penalization. Since no new data is introduced in the extrapolated region,
with probability $1 - \frac{\alpha}{2}$, where $\Phi$ denotes the cumulative distribution function of the standard normal distribution.

Figure 12 represents the results of the constrained extrapolation presented in this section, with the associated standard deviation which accounts for the *innovation error*. It is based on the data used to generate Figure 2. Without the dotted lines marking the boundaries of the initial smoothing region, it would not be possible to tell from those plots where the extrapolation starts, which is precisely the goal of this procedure. Figure 13 represents the ratio between the mortality rates obtained from the unconstrained solution of Equation 23 and the constrained solution of Equation 28. The unconstrained solution shows significant discrepancies both in the initial smoothing region and the extrapolated region compared to the constrained solution. Finally, Figure 14 represents ratios between standard deviation derived using the three presented extrapolation methods. The denominator is the constrained extrapolation method which accounts for the *innovation error*, which is used as a reference, while the numerator shows the unconstrained method as well as the constrained method which ignores the *innovation error*. The constrained extrapolation method which includes the *innovation error* has a higher standard deviation compared to the unconstrained method as it is always either equally or less smooth.

<figure><img src="image_14.png" /><figcaption>Figure 12: Constrained extrapolation of Whittaker-Henderson smoothing in the two-dimensional case. The contour lines of mortality rates, and the associated standard deviation, are depicted. The dotted lines delimitate the boundaries of the initial smoothing region.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Extrapolation of one-dimensional WH smoothing. The smoother is extrapolated on both sides of the initial observation range following a polynomial of degree q - 1 (in this case a straight line as q = 2).</figcaption></figure>

the smoothness criterion weighs more heavily in the optimization, prompting adjustments to the originally fitted values in order to produce a globally smoother estimate.

To address this, we follow the approach proposed by Carballo, Durban, and Lee (2021) and formulate a constrained optimization problem that enforces preservation of the original fitted values in the smoothing region. This is done by introducing a Lagrange multiplier $\omega$ and solving the following constrained problem:

$$ (\hat{\mathbf{y}}_+^*, \hat{\omega}) = \underset{\mathbf{\theta}_+^*, \omega}{\operatorname{argmin}} \left\{ (\mathbf{y}_+ - \mathbf{\theta}_+^*)^T W_+ (\mathbf{y}_+ - \mathbf{\theta}_+^*) + \mathbf{\theta}_+^{*T} P_+ \mathbf{\theta}_+^* + 2\omega^T (C \mathbf{\theta}_+^* - \hat{\mathbf{y}}) \right\}. $$

This optimization admits a closed-form solution for the constrained extrapolated estimator $\hat{\mathbf{y}}_+$ as a linear transformation of $\hat{\mathbf{y}}$. The derivation details are provided in Section F of the appendices. The final form is:

$$ \hat{\mathbf{y}}_+^* = Q^T \left[ \begin{matrix} I \\ -(P_+^{22})^{-1} P_+^{21} \end{matrix} \right] \hat{\mathbf{y}} $$

and the associated variance-covariance matrix is:

$$ V_+^* = Q^T \left[ \begin{matrix} V & -V P_+^{12} (P_+^{22})^{-1} \\ -(P_+^{22})^{-1} P_+^{21} V & (P_+^{22})^{-1} P_+^{21} V P_+^{12} (P_+^{22})^{-1} + (P_+^{22})^{-1} \end{matrix} \right] Q. $$
<figure><img src="image_15.png" /><figcaption>Figure 13: Ratio of mortality rates resulting from the extrapolation of the two-dimensional Whittaker-Henderson smoothing. The numerator corresponds to the unconstrained extrapolation and the denominator represents the constrained extrapolation presented in Figure 12.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Ratio of standard deviation of log-mortality rates resulting from the extrapolation of the two-dimensional Whittaker-Henderson smoothing. The numerator on the left plot corresponds to the unconstrained extrapolation and on the right plot to the constrained extrapolation that does not account for the *innovation error*. The denominator represents the constrained extrapolation presented in Figure 12.</figcaption></figure>
This formulation differs from the variance matrix of the unconstrained solution. Indeed, it enforces the constraint that the initial coefficients remain unchanged, as reflected by the presence of V (the original variance matrix) instead of $V_+^{11}$. The corresponding credible intervals are:

$$ \mathbb{E}(\mathbf{y}_+) | \mathbf{y}_+ \in [\hat{\mathbf{y}}_+^* \pm \Phi^{-1}(1 - \alpha/2) \sqrt{\sigma^2 \mathrm{diag}(V_+^*)}] $$

The following figures illustrate the impact of the constrained extrapolation procedure discussed above, using the LTC portfolio of 100,000 policyholders as a case study.

* Figure 9, left (mortality rates): this panel shows the estimated mortality rates obtained after applying the constrained extrapolation procedure to the two-dimensional WH smoothing model. The dotted lines indicate the boundaries of the original smoothing region. Visually, the transition from the smoothing region to the extrapolated area is seamless—the extrapolated surface naturally extends the smoothed mortality rates while respecting the original fitted values within the data range.

* Figure 9, right (standard deviation): this panel displays the posterior standard deviation (or credible interval width) associated with the extrapolated estimates. It reflects both the uncertainty from the original smoothing and the innovation error introduced in the extrapolated region. As expected, the standard deviation increases as we move away from the observed region, illustrating growing uncertainty about farther values.

<figure><img src="image_11.png" /><figcaption>Figure 9: Constrained extrapolation of 2D WH smoothing. The contour lines of mortality rates and the associated standard deviation are depicted. The dotted lines delimit the boundaries of the initial smoothing region.</figcaption></figure>

* Figure 10 (ratio of mortality rates): this heatmap shows the pointwise ratio between the unconstrained and constrained extrapolation of the mortality rates. A value above 1 indicates that the unconstrained version overshoots the constrained one at that location,
# 8 Discussion

In this paper, we have shown that Whittaker-Henderson smoothing, introduced in 1922, can be naturally interpreted in a Bayesian framework, which allows for the construction of credible intervals for it, provided that the observations are independent and their variances are known and used as weights. By linking it to the framework of duration models, we have demonstrated that in the construction of survival laws from experience, smoothing can be applied to the vector of raw exit rates using the observed number of exits as the weight vector. This is justified by the asymptotic properties of the maximum likelihood estimator of the raw exit rates. We have then established that the use of this asymptotic property comes at the cost of precision loss, especially when the number of observations is limited, and that a more precise iterative version of the smoothing can be proposed by directly solving penalized likelihood equations.

Furthermore, we have introduced an empirical Bayesian approach for selecting the smoothing parameter, based on maximizing a marginal likelihood function. The parameter selection introduces an iterative process that is added to the one generated by solving the penalized likelihood equations for a fixed smoothing parameter. In this case, there are two main possible approaches for parameter selection, depending on the order in which the iterations are nested: the *outer iteration* and the *performance iteration* approaches. The latter has significantly lower computation time compared to the former, at the expense of lacking theoretical convergence guarantees, which does not seem to pose a problem in practice. We have also introduced a rank reduction method based on the eigendecomposition of one-dimensional penalty matrices, which accelerates the selection of the smoothing parameter - our study on simulated data showed a 40-fold reduction in computation time - with negligible loss of precision.

Finally, we have addressed the issue of extrapolation of smoothing and shown that it requires solving a new optimization problem. Extrapolation is straightforward in the one-dimensional case. However, in the two-dimensional case constraints need to be imposed in order to preserve the values of the coefficients obtained during the initial smoothing step.

Whittaker-Henderson smoothing has been used for nearly a century by actuaries, with little change in the approach. In parallel, statistical theory on smoothing methods has undergone numerous developments, particularly in the past 30 years with the emergence of generalized additive models, covered notably by Hastie (2017) and Wood (2006). Our goal in this paper was to bridge the gap between these two perspectives. We also created an R package named WH which implements all the steps mentioned through the paper and should be straightforward to use. The results derived in the paper are directly applicable to other types of smoothing, such as the widely used P-splines smoothing by Eilers and Marx (1996). Compared to Whittaker-Henderson smoothing, P-splines naturally consider fewer parameters than observations, removing the need for the method introduced in Section 6.

While Whittaker-Henderson smoothing is applicable to both the one-dimensional and two-dimensional case, it does have limitations. Firstly, it requires regularly spaced observations. This aligns well with the format of life insurance pricing and reserving assumptions, which
while values below 1 indicate underestimation. We observe that discrepancies exist not only in the extrapolated region but also within the original data region—confirming that the unconstrained approach distorts the original estimates in order to achieve overall smoothness.

<figure><img src="image_12.png" /><figcaption>Figure 10: Ratio of mortality rates resulting from the extrapolation of 2D WH smoothing. The numerator corresponds to the unconstrained extrapolation and the denominator to the constrained extrapolation presented in Figure 9.</figcaption></figure>

* Figure 11 (ratio of standard deviations): this final figure includes two panels comparing uncertainty estimates.

* Left panel: ratio of standard deviation from the unconstrained extrapolation over that from the constrained extrapolation (including innovation error). The unconstrained version underestimate the actual uncertainty not only in the extrapolated region but also within the original data region, again reflecting the adjustments made to the original estimates in order to achieve overall smoothness.

* Right panel: ratio of standard deviation from the constrained extrapolation without innovation error over the fully constrained version with innovation error. This illustrates the contribution of the innovation error to the total uncertainty—it is substantial and should not be neglected.

# 8 Discussion

## Choosing the order of the penalization

Throughout this work, we have assumed second-order difference matrices for penalization. This choice is both standard and meaningful: from a Bayesian perspective, it corresponds to a prior
traditionally include age and duration spent in certain states that trigger indemnization. However, this can be a limit when the exits are not evenly distributed. For example, in the case of disability and long-term care risks, most exits occur in the first few months following entry into the state. In such situations, using a spline basis that is arranged to prioritize areas with more observations, as proposed by Wood (2017), could yield better results. Lastly, Whittaker-Henderson smoothing does not allow for the incorporation of additional explanatory variables in the experience table, starting with gender, which plays a major role in most biometric risks. Fortunately, it is possible to introduce these additional variables as random effects in the model by adopting a Smoothing Splines ANOVA approach, as described in Lee and Durban (2011) and Gu (2013), yielding results that are both precise and robust compared to a unisex or stratified model.

## References

Akaike, Hirotsugu. 1973. "Information Theory and an Extension of the Maximum Likelihood Principle." In *2nd International Symposium on Information Theory, 1973*.

Anderssen, RS, and Peter Bloomfield. 1974. "A Time Series Approach to Numerical Differentiation." *Technometrics* 16 (1): 69-75.

Brent, Richard P. 1973. "Algorithms for Minimization Without Derivatives, Chap. 4." Prentice-Hall, Englewood Cliffs, NJ.

Carballo, Alba, Maria Durban, Göran Kauermann, and Dae-Jin Lee. 2021. "A General Framework for Prediction in Penalized Regression." *Statistical Modelling* 21 (4): 293-312.

Carballo, Alba, Maria Durban, and Dae-Jin Lee. 2021. "Out-of-Sample Prediction in Multidimensional p-Spline Models." *Mathematics* 9 (15): 1761.

Eilers, Paul H. C., and Brian D. Marx. 1996. "Flexible Smoothing with B-Splines and Penalties." *Statistical Science* 11 (2): 89-102.

Gu, Chong. 1992. "Cross-Validating Non-Gaussian Data." *Journal of Computational and Graphical Statistics* 1 (2): 169-79.

———. 2013. *Smoothing Spline ANOVA Models*. Vol. 297. Springer.

Hastie, Trevor J. 2017. "Generalized Additive Models." In *Statistical Models in s*, 249-307. Routledge.

Henderson, Robert. 1924. "A New Method of Graduation." *Transactions of the Actuarial Society of America* 25: 29-40.

Kauermann, Göran. 2005. "A Note on Smoothing Parameter Selection for Penalized Spline Smoothing." *Journal of Statistical Planning and Inference* 127 (1-2): 53-69.

Lee, Dae-Jin, and Maria Durban. 2011. "P-Spline ANOVA-Type Interaction Models for Spatio-Temporal Smoothing." *Statistical Modelling* 11 (1): 49-69.

Marra, Giampiero, and Simon N Wood. 2012. "Coverage Properties of Confidence Intervals for Generalized Additive Model Components." *Scandinavian Journal of Statistics* 39 (1): 53-74.

Nelder, John A, and Roger Mead. 1965. "A Simplex Method for Function Minimization." *The Computer Journal* 7 (4): 308-13.
<figure><img src="image_13.png" /><figcaption>Figure 11: Ratio of standard deviation of log-mortality rates from the three extrapolation methods. Left: unconstrained vs constrained with innovation error. Right: constrained without vs with innovation error. In both, the denominator is the fully constrained method of Figure 9.</figcaption></figure>

belief that the log-transformed quantity of interest evolves linearly, which implies exponential behaviour on the original scale—consistent with actuarial models such as Gompertz.

The difference order directly shapes both the estimated trend and its extrapolation: higher-order penalties allow for more flexibility, but may induce unstable or erratic behaviour outside the data range. While Whittaker originally used third-order differences and higher orders can marginally improve model fit according to information criteria such as AIC, second-order penalties typically offer a robust compromise between smoothness, interpretability, and extrapolation stability. A detailed evaluation is provided in Section G of the appendices.

## Summary of contributions

This paper revisits the classical Whittaker-Henderson (WH) smoothing approach through the lens of modern statistical modelling. Each section brought forward a key practical insight:

* Section 2 established that WH smoothing is more than an empirical method. It has a firm Bayesian foundation. Under Gaussian assumptions, credibility intervals may be derived and used as practical substitutes to confidence intervals.

* Section 3 clarified how to construct observation and weight vectors in survival analysis models: using log-crude rates as observations and event counts as weights yields a sound statistical formulation.
Nelder, John Ashworth, and Robert WM Wedderburn. 1972. "Generalized Linear Models." *Journal of the Royal Statistical Society: Series A (General)* 135 (3): 370–84.

Patterson H. D., Thompson R. 1971. "Recovery of Inter-Block Information When Block Sizes Are Unequal." *Biometrika* 58: 545–54.

Reiss, Philip T, and R Todd Ogden. 2009. "Smoothing Parameter Selection for a Class of Semiparametric Linear Models." *Journal of the Royal Statistical Society: Series B (Statistical Methodology)* 71 (2): 505–23.

Schwarz, Gideon. 1978. "Estimating the Dimension of a Model." *The Annals of Statistics*, 461–64.

Wahba, Grace. 1980. *Spline Bases, Regularization, and Generalized Cross Validation for Solving Approximation Problems with Large Quantities of Noisy Data*. University of WISCONSIN.

———. 1985. "A Comparison of GCV and GML for Choosing the Smoothing Parameter in the Generalized Spline Smoothing Problem." *The Annals of Statistics*, 1378–1402.

Whittaker, Edmund T. 1922. "On a New Method of Graduation." *Proceedings of the Edinburgh Mathematical Society* 41: 63–75.

Wood, Simon N. 2006. *Generalized Additive Models: An Introduction with r*. chapman; hall/CRC.

———. 2011. "Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation of Semiparametric Generalized Linear Models." *Journal of the Royal Statistical Society: Series B (Statistical Methodology)* 73 (1): 3–36.

———. 2017. "P-Splines with Derivative Based Penalties and Tensor Product Smoothing of Unevenly Distributed Data." *Statistics and Computing* 27: 985–89.