A Fast Optimization Method for Additive Model via Partial Generalized Ridge Regression

Abstract

Although many statistical tools for analyzing the additive model apply the back-fitting algorithm, it is well-known that it is not guaranteed that the algorithm will converge or obtain a unique solution. Furthermore, running the algorithm on large datasets is computationally challenging. In order to address these issues, we propose a new optimization method for the additive model via partial generalized ridge regression. Using the proposed method, all trends of the additive model are estimated simultaneously, and closed-form expressions for the smoothing parameters that minimize the GCV criterion are derived. From a numerical study, compared to the back-fitting algorithm, our new method achieved higher performance with respect to both predictive accuracy and computational efficiency.

Appendix: The Proof of Equation (13)

A singular value decomposition of $(\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1}) {\varvec{W}}_2$ is expressed as

$$\begin{aligned} (\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1}) {\varvec{W}}_2 = \varvec{G} \left( \begin{array}{c} \varvec{D}^{1/2}\\ \varvec{O}_{n-k,k} \end{array}\right) \varvec{Q}' = {\varvec{G}}_1 {\varvec{D}}^{1/2} {\varvec{Q}}', \end{aligned}$$

(14)

where $\varvec{G}$ is an $n\times n$ orthogonal matrix and ${\varvec{G}}_1$ is an $n \times k$ matrix derived from the partition ${\varvec{G}}= ({\varvec{G}}_1, {\varvec{G}}_2)$. Note that ${\varvec{D}}^{1/2}$ is diagonal matrix. From Eqs. (8) and (14), we can see that

$$\begin{aligned} {\varvec{H}}_{{\varvec{\Lambda }}}&= \varvec{P}_{{\varvec{W}}_1} + {\varvec{G}}\left( \begin{array}{cc} (\varvec{D} + \varvec{\Lambda })^{-1} \varvec{D} &{}\varvec{O}_{k,n-k} \\ \varvec{O}_{n-k,k} &{}\varvec{O}_{n-k,n-k} \end{array}\right) {\varvec{G}}'. \end{aligned}$$

Hence, $\mathrm{tr}({\varvec{H}}_{{\varvec{\Lambda }}})$ can be calculated as

$$\begin{aligned} \mathrm{tr}({\varvec{H}}_{{\varvec{\Lambda }}}) = \mathrm{tr}({\varvec{P}}_{{\varvec{W}}1}) + \mathrm{tr}\{({\varvec{D}}+{\varvec{\Lambda }})^{-1} {\varvec{D}}\} = 3p + k + 1 - \sum _{j=1}^k \frac{\lambda _j}{d_j + \lambda _j}. \end{aligned}$$

(15)

Notice that ${\varvec{G}}_1 = ({\varvec{I}}_n - {\varvec{P}}_{{\varvec{W}}_1}) {\varvec{W}}_2 {\varvec{Q}}{\varvec{D}}^{-1/2}$ from (14), and $\varvec{I}_n = \varvec{G} \varvec{G}' = \varvec{G}_1 \varvec{G}_1'+\varvec{G}_2 \varvec{G}_2'$. Hence, we have $\varvec{G}_1 ' (\varvec{I}_n - {\varvec{P}}_{{\varvec{W}}_1}){\varvec{y}}= {\varvec{z}}$ and

$$\begin{aligned}&\quad ~ {\varvec{y}}'(\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1})\varvec{G}_2 \varvec{G}_2'(\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1}){\varvec{y}}\nonumber \\&= {\varvec{y}}'(\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1})(\varvec{I}_n - \varvec{G}_1 \varvec{G}_1')(\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1}){\varvec{y}}\nonumber \\&= {\varvec{y}}'\left\{ \varvec{I}_n - \varvec{P}_{{\varvec{W}}_1}- (\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1}) {\varvec{W}}_2 \varvec{Q} \varvec{D}^{-1}\varvec{Q}' {\varvec{W}}_2'(\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1})\right\} {\varvec{y}}\nonumber \\&= {\varvec{y}}' ({\varvec{I}}_n - {\varvec{H}}){\varvec{y}}= (n - 3p - k - 1) s_0^2. \end{aligned}$$

Using the above results and noting $P_{W_1}P_1=O_{n,k}$, we can derive the following equation:

$$\begin{aligned}&\quad ~ \Vert (\varvec{I}_n-\varvec{H}_{\varvec{\Lambda }}) \varvec{y}\Vert ^2 \nonumber \\&= \varvec{y} '(\varvec{I}_n - \varvec{H}_{\varvec{\Lambda }}) ^2 \varvec{y} \nonumber \\&= \varvec{y}'(\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1}) {\varvec{G}}\left\{ \varvec{I}_n - \left( \begin{array}{cc} (\varvec{D} + \varvec{\Lambda })^{-1} \varvec{D} &{}\varvec{O}_{k,n-k}\\ \varvec{O}_{n-k,k} &{}\varvec{O}_{n-k,n-k} \end{array} \right) \right\} ^2 {\varvec{G}}' (\varvec{I}_n - \varvec{P}_{{\varvec{W}}_1})\varvec{y} \nonumber \\&= {\varvec{z}}'\{ \varvec{I}_k - (\varvec{D} + \varvec{\Lambda })^{-1} \varvec{D} \}^2 {\varvec{z}}+ (n - 3p - k-1) s_0^2 \nonumber \\&= (n-3p - k-1) s_0^2+ \sum _{j=1}^k \left( \frac{\lambda _j}{\lambda _j + d_j}z_j\right) ^2. \end{aligned}$$

By substituting (15) and above result into (9), we can obtain

$$\begin{aligned} \mathrm{GCV}({\varvec{\Lambda }}) = \dfrac{(n - 3p - k - 1) s_0^2 + \sum _{j=1}^{k} \{ \lambda _j z_j / (\lambda _j + d_j) \}^2 }{n [1 - \{ 3p + k + 1 -\sum _{j=1}^k \lambda _j / (\lambda _j + d_j) \}/n ] ^2 } . \end{aligned}$$

(16)

Let ${\varvec{\delta }}= (\delta _1, \dots , \delta _k)'$ be a k-dimensional vector defined as

$$ \delta _j = \frac{\lambda _j}{d_j + \lambda _j}. $$

Note that $\delta _j \in [0,1]$ since $d_j \ge 0$ and $\lambda _j \ge 0$. Then, the GCV criterion in (16) is expressed as the following function with respect to ${\varvec{\delta }}$:

$$\begin{aligned} \mathrm{GCV}({\varvec{\Lambda }}) = f ({\varvec{\delta }}) = \frac{r({\varvec{\delta }})}{c({\varvec{\delta }})^2}, \end{aligned}$$

(17)

where the functions $r({\varvec{\delta }})$ and $c({\varvec{\delta }})$ are given by

$$\begin{aligned} r({\varvec{\delta }}) = \frac{(n-3p - k-1) s_0^2 + \sum _{j=1}^k z_j^2 \delta _j^2 }{n}, \ c({\varvec{\delta }}) = 1- \frac{1}{n}\left( 3p + k + 1 - \sum _{j=1}^k \delta _j\right) . \end{aligned}$$

Here, $z_1,\dots ,z_k$ are given in (10). Let $\hat{{\varvec{\delta }}} = (\hat{\delta }_1, \dots , \hat{\delta }_k)'$ be the minimizer of $f({\varvec{\delta }})$ in (17), i.e.,

$$ \hat{{\varvec{\delta }}} = \arg \min _{{\varvec{\delta }}\in [0,1]^k} f({\varvec{\delta }}), $$

where $[0,1]^k$ is the kth Cartesian power of the set [0, 1]. Notice that $r(\delta )$ and $c(\delta )$ are differentiable functions with respect to $\delta _j$. Thus, we obtain

$$ \frac{\partial }{\partial \delta _j} f({\varvec{\delta }}) = \frac{2}{nc({\varvec{\delta }})^3}\left\{ c({\varvec{\delta }}) z_j^2 \delta _j - r({\varvec{\delta }})\right\} . $$

Hence, noting $c(\delta )$ is a finite function, we find a necessary condition for $\hat{{\varvec{\delta }}}$ as

$$\begin{aligned} \hat{\delta }_j = \left\{ \begin{array}{ll} 1 &{} (h(\hat{{\varvec{\delta }}}) \ge z_j^2) \\ h(\hat{{\varvec{\delta }}})/z_j^2 &{}(h(\hat{{\varvec{\delta }}}) < z_j^2) \end{array} \right. , \end{aligned}$$

where $h({\varvec{\delta }}) = r({\varvec{\delta }})/c({\varvec{\delta }}) > 0$.

On the other hand, let ${\varvec{H}}= \{ {\varvec{\delta }}\in [0, 1]^k| {\varvec{\delta }}= {\varvec{\delta }}^{\star } (h), \forall h \in \mathbb {R}_+ \}$, where ${\varvec{\delta }}^\star (h)$ is the k-dimensional vector for which the jth element is defined as

$$ \delta _j^\star (h) = \left\{ \begin{array}{ll} 1 &{} (h \ge z_j^2) \\ h /z_j^2 &{}(h < z_j^2) \end{array} \right. , $$

and $\mathbb {R}_+$ is a set of nonnegative real numbers. Then, it follows from $\mathcal {H}\subseteq [0, 1]^k$ and $\hat{{\varvec{\delta }}} \in \mathcal {H}$ that

$$ f (\hat{{\varvec{\delta }}}) = \min _{{\varvec{\delta }}\in [0, 1]^k} f ({\varvec{\delta }}) \le \min _{{\varvec{\delta }}\in \mathcal {H}} f ({\varvec{\delta }}) = \min _{h \in \mathbb {R}_+} f ({\varvec{\delta }}^\star (h)), \quad f (\hat{{\varvec{\delta }}}) \ge \min _{h \in \mathbb {R}_+} f ({\varvec{\delta }}^\star (h)). $$

Hence, we have

$$\begin{aligned} \hat{{\varvec{\delta }}} = {\varvec{\delta }}^\star (\hat{h})\quad \left( \hat{h} = \arg \min _{h \in \mathbb {R}_+} f ({\varvec{\delta }}^\star (h)) \right) . \end{aligned}$$

Using this result, we minimize the following function:

$$\begin{aligned} f ({\varvec{\delta }}^\star (h)) = f_1 (h) = \dfrac{ r_1 (h) }{ c_1 (h)^2 }, \end{aligned}$$

where $r_1(h)=r({\varvec{\delta }}^\star (h))$ and $c_1(h)=c({\varvec{\delta }}^\star (h))$, which can be calculated as

$$\begin{aligned} r_1 (h)&= \frac{1}{n} \left[ (n-3p - k-1) s_0^2 + \sum _{j=1}^k \left\{ I (h<z_j^2 ) \left( \dfrac{h}{z_j^2} - 1 \right) + 1 \right\} ^2 z_j^2 \right] , \\ c_1 (h)&= 1- \frac{1}{n} \left[ 3p + k + 1 - \sum _{j=1}^k \left\{ I (h <z_j^2 ) \left( \dfrac{h}{z_j^2} - 1 \right) + 1 \right\} \right] . \end{aligned}$$

Suppose that $h \in R_a$, where $a \in \{0,1,\ldots ,k\}$ and $R_a$ is a range defined by (12). Notice that $u_j$ are the jth-order statistics of $z_1^2, \ldots , z_k^2$. Then, we have

$$ f_1 (h) = f_{1a} (h) = \dfrac{ r_{1a} (h) }{ c_{1a} (h)^2 }\quad (h \in R_a), $$

where functions $r_{1a}(h)$ and $c_{1a}(h)$ are given by

$$\begin{aligned} r_1 (h)&= r_{1a} (h) = \frac{1}{n} \left\{ (n-3p - k-1+a) s_a^2 + \ell _a h^2 \right\} , \\ c_1 (h)&= c_{1a} (h) = 1 - \frac{1}{n} ( 3p + k + 1 - a - \ell _a h ), \end{aligned}$$

where $s_a^2$ is given in (11) and $\ell _a = \sum _{j=a+1}^k 1/u_j$. By simple calculation, let $g_a (h) = (n - 3p - k - 1 + a) (h - s_a^2)$, we can obtain

$$\begin{aligned} \dfrac{d}{dh} f_{1a} (h)&= \dfrac{2 \ell _a}{n^2 c_{1a} (h)^3} g_a (h). \end{aligned}$$

Here, we note $g_a (u_{a+1})= g_{a+1} (u_{a+1}) \quad (a = 0, \dots , k-1).$ Moreover, $2 \ell _a/\{n^2 c_{1a} (h)^3\}$ is positive, $\lim _{h \rightarrow 0} g_0 (h) < 0$, and $g_a (h)\ (h \in R_a)$ is a monotonically increasing function in $h \in \mathbb {R}_+$, since $n-3p-k-1>0$ and $a \ge 0$. Consequently, the Eq. (13) is obtained by combining $\hat{\delta }_j$ and $h^*=s_a^2$ where $d f_{1a}(h)/d h|_{h=h^*}=0$, and some calculation.