simulates the overall proportion of heads in a run of $n$ separate coin flips.

######  Definition 8.1 (Binomial Distribution)

A random variable $X$ is said to follow the binomial distribution with
parameter $p\in(0,1)$ and $n\in\mathbb{N}$, denoted by
$X\sim\mathrm{Binom}(n,p)$ if

$f(x;n,p)=\binom{n}{x}p^{x}(1-p)^{n-x}.$

The mean and variance of $X\sim\mathrm{Binom}(a,b)$ are given by

$\mathrm{E}[X]=np,\qquad\mathrm{Var}[X]=np(1-p).$

Figure 8.1 compares different parameters of $p$ with $n=10$ for the binomial
distribution.

Figure 8.1: Binomial distribution probability mass functions for different
values of the parameters $p$ with $n=10$.

And we will see the prior under this model is the probability density function
of Beta distribution:

$\mathrm{prior}=\mathrm{Beta}(\theta|a,b)=p(\theta|a,b)=\frac{1}{B(a,b)}\theta^{a-1}(1-\theta)^{b-1}\mathds{1}(0\leq\theta\leq
1),$

where $B(a,b)$ is the Euler’s beta function and it can be seen as a
normalization term.

######  Definition 8.2 (Beta Distribution)

A random variable $X$ is said to follow the beta distribution with parameter
$a>0$ and $b>0$, denoted by $X\sim\mathrm{Beta}(a,b)$ if

$f(x;a,b)=\left\\{\begin{aligned}
&\frac{1}{B(a,b)}x^{a-1}(1-x)^{b-1},&\mathrm{\,\,if\,\,}0\leq x\leq 1.\\\
&0,&\mathrm{\,\,otherwise\,\,},\end{aligned}\right.$

where $B(a,b)$ is the Euler’s beta function and it can be seen as a
normalization term. Equivalently, $B(a,b)$ can be obtained by

$B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)},$

where $\Gamma(\cdot)$ is the gamma function (Definition 3.1, p. 3.1). The mean
and variance of $X\sim\mathrm{Beta}(a,b)$ are given by

$\mathrm{E}[X]=\frac{a}{a+b},\qquad\mathrm{Var}[X]=\frac{ab}{(a+b+1)(a+b)^{2}}.$

Figure 8.2 compares different parameters of $a,b$ for beta distribution. When
$a=b=1$, the beta distribution is just a uniform distribution in the range of
0 and 1.

We put a Beta prior over the parameter $\theta$ of the Bernoulli distribution.
The posterior is obtained by

$\displaystyle\mathrm{posterior}=p(\theta|x_{1:n})$ $\displaystyle\propto
p(x_{1:n}|\theta)p(\theta|a,b)$ $\displaystyle=\theta^{\sum
x_{i}}(1-\theta)^{n-\sum
x_{i}}\times\frac{1}{B(a,b)}\theta^{a-1}(1-\theta)^{b-1}\mathds{1}(0<\theta<1)$
$\displaystyle\propto\theta^{a+\sum x_{i}-1}(1-\theta)^{b+n-\sum
x_{i}-1}\mathds{1}(0<\theta<1)$ $\displaystyle\propto Beta(\theta|a+\sum
x_{i},b+n-\sum x_{i}).$

Figure 8.2: Beta distribution probability density functions for different
values of the parameters $a,b$. When $a=b=1$, the beta distribution is the
uniform distribution in the range of 0 and 1.

We find that the posterior distribution shares the same form as the prior
distribution. When this happens, we call the prior conjugate prior. The
conjugate prior has a nice form such that it is easy to work with for
computing the posterior probability density function and its derivatives,
sampling from the posterior.

######  Remark 8.1 (Prior Information in Beta-Bernoulli Model)

A comparison of the prior and posterior formulation would find that the
hyperparameter $a$ is the prior number of $1$’s in the output and $b$ is the
prior number of 0’s in the output. And $a+b$ is the prior information about
the sample size.

######  Remark 8.2 (Bayesian Estimator)

From this example of the Beta-Bernoulli model, like maximum likelihood
estimator and method of moment (MoM, i.e., using the moment information to get
the model parameter.), the Bayesian model is also a kind of point estimator.
But Bayesian models output a probability of the parameter of interest
$p(\theta|x_{1:n})$.

When we want to predict for new coming data, we do not give out the prediction
by a direct model $p(x_{n+1}|\theta)$. But rather an integration:

$p(x_{n+1}|x_{1:n})=\int p(x_{n+1}|\theta)p(\theta|x_{1:n})d\theta.$

In another word, $x_{n+1}$ is dependent on $x_{1:n}$. $x_{1:n}$ provide
information on $\theta$, which in turn provides information on $x_{n+1}$
(i.e., $x_{1:n}\rightarrow\theta\rightarrow x_{n+1}$).

######  Example 8.2.1 (Amount of Data Matters)

Suppose we have three observations for the success in Bernoulli distribution:

1). 10 out of 10 are observed to be success (1’s);

2). 48 out of 50 are observed to be success (1’s);

3). 186 out of 200 are observed to be success (1’s).

So, what is the probability of success in the Bernoulli model? Normal answer
to case 1, 2, 3 are 100%, 96% and 93% respectively. But an observation of 10
inputs is rather a small amount of data and noise can make it less convincing.

Suppose we put a $\mathrm{Beta}(1,1)$ (a uniform distribution, see Figure 8.2)
prior over the Bernoulli distribution. The posterior probability of success
for each case would be $\frac{11}{12}=91.6\%$, $\frac{49}{52}=94.2\%$ and
$\frac{187}{202}=92.6\%$ respectively. Now we find the case 1 has less
probability of success compared to case 2.

A Bayesian view of the problem naturally incorporates the amount of data as
well as its average. This special case shown here is also called the Laplace’s
rate of succession (Ollivier, 2015). Laplace’s “add-one” rule of succession
modifies the observed frequencies in a sequence of successes and failures by
adding one to the observed counts. This improves prediction by avoiding zero
probabilities and corresponds to a uniform Bayesian prior on the parameter.
$\square$

This example above shows that Bayesian models consider prior information on
the parameters in the model making it particularly useful to regularize
regression problems where data information is limited. And this is why the
Bayesian approach gains worldwide attention for decades. The prior information
$p(\theta)$ and likelihood function $p(x|\theta)$ represent a rational
person’s belief, and then the Bayes’ rule is an optimal method of updating
this person’s beliefs about $\theta$ given new information from the data
(Fahrmeir et al., 2007; Hoff, 2009). The prior information given by
$p(\theta)$ might be wrong if it does not accurately represent our prior
beliefs. However, this does not mean that the posterior $p(\theta|x)$ is not
useful. A famous quote is “all models are wrong, but some are useful” (Box and
Draper, 1987). If the prior $p(\theta)$ approximates our beliefs, then the
posterior $p(\theta|x)$ is also a good approximation to posterior beliefs.

###  8.3 Bayesian Linear Model: Zero-Mean Prior

We now present the Bayesian methods for linear models. Assume
$\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}$ where
$\boldsymbol{\epsilon}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\boldsymbol{I})$
and $\sigma^{2}$ is fixed. As discussed in Section 2 (p. 2), this additive
Gaussian noise assumption gives rise to the likelihood. Let
$\mathcal{X}(\bm{x}_{1:n})=\\{\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{n}\\}$ be
the observations of $n$ data points,

$\mathrm{likelihood}=\bm{y}|\bm{X},\boldsymbol{\beta},\sigma^{2}\sim\mathcal{N}(\bm{X}\boldsymbol{\beta},\sigma^{2}\bm{I}).$

Suppose we specify a Gaussian prior with zero-mean over the weight parameter

$\mathrm{prior}=\boldsymbol{\beta}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{0}).$

By the Bayes’ theorem
“$\mathrm{posterior}\propto\mathrm{likelihood}\times\mathrm{prior}$”, we get
the posterior

$\displaystyle\mathrm{posterior}$
$\displaystyle=p(\boldsymbol{\beta}|\bm{y},\bm{X},\sigma^{2})$
$\displaystyle\propto
p(\bm{y}|\bm{X},\boldsymbol{\beta},\sigma^{2})p(\boldsymbol{\beta}|\boldsymbol{\Sigma}_{0})$
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{n/2}}\exp\left(-\frac{1}{2\sigma^{2}}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})\right)$
$\displaystyle\,\,\,\,\,\,\,\times\frac{1}{(2\pi)^{n/2}|\boldsymbol{\Sigma}_{0}|^{1/2}}\exp\left(-\frac{1}{2}\boldsymbol{\beta}^{\top}\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}\right)$
$\displaystyle\propto\exp\left(-\frac{1}{2}(\boldsymbol{\beta}-\boldsymbol{\beta}_{1})^{\top}\boldsymbol{\Sigma}_{1}^{-1}(\boldsymbol{\beta}-\boldsymbol{\beta}_{1})\right),$

where
$\boldsymbol{\Sigma}_{1}=(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}$
and
$\boldsymbol{\beta}_{1}=(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{y})$.
Therefore the posterior distribution is also a Gaussian distribution (same
form as the prior distribution):

$\mathrm{posterior}=\boldsymbol{\beta}|\bm{y},\bm{X},\sigma^{2}\sim\mathcal{N}(\boldsymbol{\beta}_{1},\boldsymbol{\Sigma}_{1}).$

##### A word on the notation

Note that we use $\\{\boldsymbol{\beta}_{1},\boldsymbol{\Sigma}_{1}\\}$ to
denote the posterior mean and posterior covariance in the zero-mean prior
model. Similarly, the posterior mean and posterior covariance in semi-
conjugate prior and full-conjugate prior models will be denoted as
$\\{\boldsymbol{\beta}_{2},\boldsymbol{\Sigma}_{2}\\}$ and
$\\{\boldsymbol{\beta}_{3},\boldsymbol{\Sigma}_{3}\\}$ respectively (see
sections below).

In this case, we do not need to assume $\boldsymbol{X}$ has full rank
generally. Note further that if we assume $\bm{X}$ has full rank, in the
limit, when $\boldsymbol{\Sigma}_{0}\rightarrow\boldsymbol{0}$,
$\boldsymbol{\beta}_{1}\rightarrow\hat{\boldsymbol{\beta}}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}\bm{y}$,
in which case, maximum a posteriori (MAP) estimator from Bayesian model goes
back to ordinary least squares estimator. And the posterior is
$\boldsymbol{\beta}|\bm{y},\bm{X},\sigma^{2}\sim\mathcal{N}(\hat{\boldsymbol{\beta}},\sigma^{2}(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1})$,
which shares similar form as the OLS estimator
$\hat{\boldsymbol{\beta}}\sim\mathcal{N}(\boldsymbol{\beta},\sigma^{2}(\boldsymbol{X}^{\top}\boldsymbol{X})^{-1})$
under Gaussian disturbance (see Lemma 3.2, p. 3.2).

######  Remark 8.3 (Ridge Regression)

In least squares approximation, we use $\bm{X}\boldsymbol{\beta}$ to
approximate $\bm{y}$. Two issues arise: the model can potentially overfit and
$\bm{X}$ may not have full rank. In ridge regression, we regularize large
value of $\boldsymbol{\beta}$ and thus favor simpler models. Instead of
minimizing $||\bm{y}-\bm{X}\boldsymbol{\beta}||^{2}$, we minimize
$||\bm{y}-\bm{X}\boldsymbol{\beta}||^{2}+\lambda||\boldsymbol{\beta}||^{2}$,
where $\lambda$ is a hyper-parameter that can be tuned:

$\mathop{\arg\min}_{\boldsymbol{\beta}}{(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})+\lambda\boldsymbol{\beta}^{\top}\boldsymbol{\beta}}.$

By differentiating and setting the derivative to zero we get

$\hat{\boldsymbol{\beta}}_{ridge}=(\bm{X}^{\top}\bm{X}+\lambda\bm{I})^{-1}\bm{X}^{\top}\bm{y},$

in which case, $(\bm{X}^{\top}\bm{X}+\lambda\bm{I})$ is invertible even when
$\bm{X}$ does not have full rank. We leave more details about ridge regression
to the readers.

Realize that when we set $\boldsymbol{\Sigma}_{0}=\bm{I}$, we obtain
$\boldsymbol{\beta}_{1}=(\bm{X}^{\top}\bm{X}+\sigma^{2}\bm{I})^{-1}\bm{X}^{\top}\bm{y}$
and
$\boldsymbol{\Sigma}_{1}=(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\bm{I})^{-1}$.
Since
$\mathrm{posterior}=\boldsymbol{\beta}|\bm{y},\bm{X},\sigma^{2}\sim\mathcal{N}(\boldsymbol{\beta}_{1},\boldsymbol{\Sigma}_{1})$.
The MAP estimator of
$\boldsymbol{\beta}=\boldsymbol{\beta}_{1}=(\bm{X}^{\top}\bm{X}+\sigma^{2}\bm{I})^{-1}\bm{X}^{\top}\bm{y}$,
which shares the same form as ridge regression by letting
$\sigma^{2}=\lambda$. Thus we notice ridge regression is a special case of
Bayesian linear model with zero-mean prior. And ridge regression has a nice
interpretation from the Bayesian approach - finding the mode of the posterior.
An example is shown in Rasmussen (2003) where the “well determined” (i.e., the
distribution around the slope is more compact) slope of $\boldsymbol{\beta}$
is almost unchanged after the posterior process while the intercept (which is
more dispersed) shrunk towards zero. This is actually a regularization effect
on the parameter like ridge regression.

####  8.3.1 Zeller’s $g$-Prior and Variable Transformation

As a specific example, suppose we want to analyze a person’s weight by human’s
characteristics and a variable in $\bm{X}$ denotes the height of a person
which is in meter. What if now the variable is in centimeter? The model is the
same if we divide the related parameter in $\boldsymbol{\beta}$ by 100
(transfer centimeter back into meter measurement). More generally, suppose
input matrix $\bm{X}$ is transferred by $\tilde{\bm{X}}=\bm{X}\bm{P}$ given
some $p\times p$ matrix $\bm{P}$ in which case the model parameter is
$\tilde{\boldsymbol{\beta}}$. Then, we have

$\displaystyle\bm{y}=\bm{X}\boldsymbol{\beta}=\tilde{\bm{X}}\tilde{\boldsymbol{\beta}}=\bm{X}\bm{P}\tilde{\boldsymbol{\beta}}.$

According to the principle of invariance, the posterior distributions of
$\boldsymbol{\beta}$ and $\bm{P}\tilde{\boldsymbol{\beta}}$ should be the
same. As shown previously, for $\bm{X}$, the posterior is

$\boldsymbol{\beta}|\bm{y},\bm{X},\sigma^{2}\sim\mathcal{N}\left((\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{y}),(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}\right).$

Similarly, for $\tilde{\bm{X}}$, the posterior
$\bm{P}\tilde{\boldsymbol{\beta}}$ is

$\bm{P}\tilde{\boldsymbol{\beta}}|\bm{y},\tilde{\bm{X}},\sigma^{2}\sim\mathcal{N}\left({\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{P}}(\frac{1}{\sigma^{2}}\tilde{\bm{X}}^{\top}\tilde{\bm{X}}+\tilde{\boldsymbol{\Sigma}_{0}}^{-1})^{-1}(\frac{1}{\sigma^{2}}\tilde{\bm{X}}^{\top}\bm{y}),{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{P}}(\frac{1}{\sigma^{2}}\tilde{\bm{X}}^{\top}\tilde{\bm{X}}+\tilde{\boldsymbol{\Sigma}_{0}}^{-1})^{-1}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{P}^{\top}}\right).$

A simple calculation can show that this condition is met if
$\boldsymbol{\Sigma}_{0}=k(\bm{X}^{\top}\bm{X})^{-1}$ where $k>0$ is a
hyperparameter. A popular specification of $k$ is to relate it to the noise
variance $\sigma^{2}$ by $k=g\sigma^{2}$. This is called the Zeller’s
$g$-prior (Zellner, 1986). Following the Bayesian linear model with zero-mean
prior, the posterior of $\boldsymbol{\beta}$ is

$\mathrm{posterior}=\boldsymbol{\beta}|\bm{y},\bm{X},\sigma^{2}\sim\mathcal{N}(\boldsymbol{\beta}_{1},\boldsymbol{\Sigma}_{1}).$

where
$\boldsymbol{\Sigma}_{1}=\frac{g\sigma^{2}}{g+1}(\bm{X}^{\top}\bm{X})^{-1}$
and
$\boldsymbol{\beta}_{1}=\frac{g}{g+1}(\bm{X}^{\top}\bm{X})^{-1}(\bm{X}^{\top}\bm{y})$.

###  8.4 Bayesian Linear Model: Semiconjugate Prior Distribution

We will use gamma distribution as the prior of the inverse variance
(precision) parameter in Gaussian distribution.

######  Definition 8.3 (Gamma Distribution)

A random variable $X$ is said to follow the gamma distribution with parameter
$r>0$ and $\lambda>0$, denoted by $X\sim\mathrm{Gamma}(r,\lambda)$ if

$f(x;r,\lambda)=\left\\{\begin{aligned}
&\frac{\lambda^{r}}{\Gamma(r)}x^{r-1}exp(-\lambda x),&\mathrm{\,\,if\,\,}x\geq
0.\\\ &0,&\mathrm{\,\,if\,\,}x<0.\end{aligned}\right.$

So if $X\sim\chi_{(p)}^{2}$, then $X\sim\mathrm{Gamma}(p/2,1/2)$, i.e., Chi-
square distribution is a special case of the gamma distribution. The mean and
variance of $X\sim\mathrm{Gamma}(r,\lambda)$ are given by

$\mathrm{E}[X]=\frac{r}{\lambda},\qquad\mathrm{Var}[X]=\frac{r}{\lambda^{2}}.$

Figure 8.3 compares different parameters for gamma distribution and Chi-square
distribution.

(a) Gamma probability density functions for different values of the parameters
$r,\lambda$.

(b) Chi-square probability density functions for different values of the
parameter $p$.

Figure 8.3: Comparison between the gamma distribution and Chi-square
distribution (Definition 3.1, p. 3.1).

As for the reason of using the gamma distribution as the prior for precision,
we quote the description from Kruschke (2014):  Because of its role in
conjugate priors for normal likelihood function, the gamma distribution is
routinely used as a prior for precision (i.e., inverse variance). But there is
no logical necessity to do so, and modern MCMC methods permit more flexible
specification of priors. Indeed, because precision is less intuitive than
standard deviation, it can be more useful to give standard deviation a uniform
prior that spans a wide range.

Same setting as Section 8.3 (p. 8.3), but we assume now $\sigma^{2}$ is not
fixed. Again, we have the likelihood function by

$\mathrm{likelihood}=\bm{y}|\bm{X},\boldsymbol{\beta},\sigma^{2}\sim\mathcal{N}(\bm{X}\boldsymbol{\beta},\sigma^{2}\bm{I}).$

We specify a non zero-mean Gaussian prior over the weight parameter

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathrm{prior:\,}}$
$\displaystyle\boldsymbol{\beta}\sim\mathcal{N}({\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{\beta}_{0}},\boldsymbol{\Sigma}_{0})$
$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\gamma=1/\sigma^{2}\sim\mathrm{Gamma}(a_{0},b_{0})},$

where we differentiate from previous descriptions by blue text.

(1). Then, given $\sigma^{2}$, by the Bayes’ theorem
“$\mathrm{posterior}\propto\mathrm{likelihood}\times\mathrm{prior}$”, we get
the posterior

$\displaystyle\mathrm{posterior}$
$\displaystyle=p(\boldsymbol{\beta}|\bm{y},\bm{X},\sigma^{2})$
$\displaystyle\propto
p(\bm{y}|\bm{X},\boldsymbol{\beta},\sigma^{2})p(\boldsymbol{\beta}|\boldsymbol{\beta}_{0},\boldsymbol{\Sigma}_{0})$
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{n/2}}\exp\left(-\frac{1}{2\sigma^{2}}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})\right)$
$\displaystyle\,\,\,\,\,\,\times\frac{1}{(2\pi)^{n/2}|\boldsymbol{\Sigma}_{0}|^{1/2}}\exp\left(-\frac{1}{2}(\boldsymbol{\beta}-\boldsymbol{\beta}_{0})^{\top}\boldsymbol{\Sigma}_{0}^{-1}(\boldsymbol{\beta}-\boldsymbol{\beta}_{0})\right)$
$\displaystyle\propto\exp\left(-\frac{1}{2}(\boldsymbol{\beta}-\boldsymbol{\beta}_{2})^{\top}\boldsymbol{\Sigma}_{2}^{-1}(\boldsymbol{\beta}-\boldsymbol{\beta}_{2})\right),$

where
$\boldsymbol{\Sigma}_{2}=(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}$
and

$\boldsymbol{\beta}_{2}=\boldsymbol{\Sigma}_{2}(\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0}+\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{y})=(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}({\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0}}+\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{y}).$

Therefore, the posterior density is also a Gaussian distribution:

$\mathrm{posterior}=\boldsymbol{\beta}|\bm{y},\bm{X},\sigma^{2}\sim\mathcal{N}(\boldsymbol{\beta}_{2},\boldsymbol{\Sigma}_{2}).$

1\. $\boldsymbol{\Sigma}_{0}$ here is a fixed hyperparameter. 2\. We note that
$\boldsymbol{\beta}_{1}$ in Section 8.3 (p. 8.3) is a special case of
$\boldsymbol{\beta}_{2}$ when $\boldsymbol{\beta}_{0}=\boldsymbol{0}$. 3\. And
if we assume further $\bm{X}$ has full rank. When
$\boldsymbol{\Sigma}_{0}^{-1}\rightarrow\boldsymbol{0}$,
$\boldsymbol{\beta}_{2}\rightarrow\hat{\boldsymbol{\beta}}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}\bm{y}$
which is the OLS estimator. 4\. When $\sigma^{2}\rightarrow\infty$,
$\boldsymbol{\beta}_{2}$ is approximately approaching to
$\boldsymbol{\beta}_{0}$, the prior expectation of parameter. However, in
zero-mean prior, $\sigma^{2}\rightarrow\infty$ will make
$\boldsymbol{\beta}_{1}$ approach to $\boldsymbol{0}$. 5\. Weighted average:
we reformulate by $\displaystyle\boldsymbol{\beta}_{2}$
$\displaystyle=(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}(\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0}+\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{y})$
$\displaystyle=(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0}+(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}\frac{\bm{X}^{\top}\bm{X}}{\sigma^{2}}(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{y}$
$\displaystyle=(\bm{I}-\bm{A})\boldsymbol{\beta}_{0}+\bm{A}\hat{\boldsymbol{\beta}},$
where $\hat{\boldsymbol{\beta}}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{y}$
is the OLS estimator of $\boldsymbol{\beta}$ and
$\bm{A}=(\frac{1}{\sigma^{2}}\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}\frac{\bm{X}^{\top}\bm{X}}{\sigma^{2}}$.
We see that the posterior mean of $\boldsymbol{\beta}$ is a weighted average
of prior mean and OLS estimator of $\boldsymbol{\beta}$. Thus, if we set the
prior parameter $\boldsymbol{\beta}_{0}=\hat{\boldsymbol{\beta}}$, the
posterior mean of $\boldsymbol{\beta}$ will be exactly
$\hat{\boldsymbol{\beta}}$.

(2). Given $\boldsymbol{\beta}$, again, by Bayes’ theorem, we obtain the
posterior

$\displaystyle\mathrm{posterior}$
$\displaystyle=p(\gamma=\frac{1}{\sigma^{2}}|\bm{y},\bm{X},\boldsymbol{\beta})$
$\displaystyle\propto
p(\bm{y}|\bm{X},\boldsymbol{\beta},\gamma)p(\gamma|a_{0},b_{0})$
$\displaystyle=\frac{\gamma^{n/2}}{(2\pi)^{n/2}}\exp\left(-\frac{\gamma}{2}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})\right)$
$\displaystyle\,\,\,\,\,\,\,\times\frac{{b_{0}}^{a_{0}}}{\Gamma(a_{0})}\gamma^{a_{0}-1}\exp(-b_{0}\gamma)$
$\displaystyle\propto\gamma(a_{0}+\frac{n}{2}-1)\exp\left(-\gamma\left[b_{0}+\frac{1}{2}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})\right]\right),$

and the posterior is a Gamma distribution:

$\mathrm{posterior\,\,of\,\,}\gamma\mathrm{\,\,given\,\,}\boldsymbol{\beta}=\gamma|\bm{y},\bm{X},\boldsymbol{\beta}\sim\mathrm{Gamma}\left(a_{0}+\frac{n}{2},[b_{0}+\frac{1}{2}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})]\right).$
(8.1)

1\. We notice that the prior mean and posterior mean of $\gamma$ are
$\mathrm{E}[\gamma]=\frac{a_{0}}{b_{0}}$ and
$\mathrm{E}[\gamma|\boldsymbol{\beta}]=\frac{a_{0}+\frac{n}{2}}{b_{0}+\frac{1}{2}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})}$
respectively. So the internal meaning of $2a_{0}$ is the prior sample size for
the noise variance parameter $\sigma^{2}=\frac{1}{\gamma}$. 2\. As we assume
$\bm{y}=\bm{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}$ where
$\boldsymbol{\epsilon}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\bm{I})$, then
$\frac{(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})}{\sigma^{2}}\sim\chi_{(n)}^{2}$
and
$\mathrm{E}[\frac{1}{2}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})]=\frac{n}{2}\sigma^{2}$.
So the internal meaning of $\frac{2b_{0}}{a_{0}}$ is the prior variance of the
noise. 3\. Some textbooks would write
$\gamma\sim\mathrm{Gamma}(n_{0}/2,n_{0}\sigma_{0}^{2}/2)$ to make this
explicit (in which case, $n_{0}$ is the prior sample size, and
$\sigma_{0}^{2}$ is the prior variance). But a prior in this form seems coming
from nowhere at a first glance.

####  8.4.1 Gibbs Sampling with Two Variables

Gibbs sampling was introduced by Turchin (Turchin, 1971), and later introduced
by brothers Geman in the context of image restoration (Geman and Geman, 1984).
The Geman brothers named the algorithm after the physicist J. W. Gibbs, some
eight decades after his death, in reference to an analogy between the sampling
algorithm and statistical physics.

Gibbs sampling is applicable when the joint distribution is not known
explicitly or is difficult to sample from directly, but the conditional
distribution of each variable is known and easy to sample from. A Gibbs
sampler generates a draw from the distribution of each parameter or variable
in turn, conditional on the current values of the other parameters or
variables. Therefore, a Gibbs sampler is a componentwise algorithm. In our
example from Section 8.1, given some data $\mathcal{X}$ and a probability
distribution $p(\boldsymbol{\beta}|\mathcal{X},\boldsymbol{\alpha})$
parameterized by
$\boldsymbol{\beta}=\\{\beta_{1},\beta_{2},\ldots,\beta_{p}\\}$. We can
successively draw samples from the distribution by sampling from

$\beta_{i}^{(t)}\sim
p(\beta_{i}|\boldsymbol{\beta}_{-i}^{(t-1)},\mathcal{X},\boldsymbol{\alpha}),$
(8.2)

where $\boldsymbol{\beta}_{-i}^{(t-1)}$ is all current values of
$\boldsymbol{\beta}$ in $(t-1)^{th}$ iteration except for $\beta_{i}$. If we
sample long enough, these $\beta_{i}$ values will be random samples from the
distribution $p$.

In deriving a Gibbs sampler, it is often helpful to observe that

$p(\beta_{i}\,|\,\boldsymbol{\beta}_{-i},\mathcal{X})=\frac{p(\beta_{1},\beta_{2},\ldots,\beta_{p},\mathcal{X})}{p(\boldsymbol{\beta}_{-i},\mathcal{X})}\propto
p(\beta_{1},\beta_{2},\ldots,\beta_{p},\mathcal{X}).$ (8.3)

The conditional distribution is proportional to the joint distribution. We
will get a lot of benefits from this simple observation by dropping constant
terms from the joint distribution (relative to the parameters we are
conditioned on).

Shortly, as a simplified example, given a joint probability distribution
$p(\beta_{1},\beta_{2}|\mathcal{X})$, a Gibbs sampler would draw
$p(\beta_{1}|\beta_{2},\mathcal{X})$ , then
$p(\beta_{2}|\beta_{1},\mathcal{X})$ iteratively. The procedure defines a
sequence of realization of random variables $\beta_{1}$ and $\beta_{2}$

$(\beta_{1}^{0},\beta_{2}^{0}),(\beta_{1}^{1},\beta_{2}^{1}),(\beta_{1}^{2},\beta_{2}^{2}),...$

which converges to the joint distribution $p(\beta_{1},\beta_{2})$. More
details about Gibbs sampling can be found in Turchin (1971); Geman and Geman
(1984); Müller and Quintana (2004); Rencher and Schaalje (2008); Hoff (2009);
Gelman et al. (2013); Kruschke and Liddell (2018).

By this Gibbs sampling method, we can construct a Gibbs sampler for the
Bayesian linear model with semiconjugate prior in Section 8.4:

0\. Set initial values to $\boldsymbol{\beta}$ and
$\gamma=\frac{1}{\sigma^{2}}$;

1\. update $\boldsymbol{\beta}$:
$\mathrm{posterior}=\boldsymbol{\beta}|\bm{y},\bm{X},\gamma\sim\mathcal{N}(\boldsymbol{\beta}_{2},\boldsymbol{\Sigma}_{2})$;

2\. update $\gamma$:
$\mathrm{posterior}=\gamma|\bm{y},\bm{X},\boldsymbol{\beta}\sim\mathrm{Gamma}\left(a_{0}+\frac{n}{2},[b_{0}+\frac{1}{2}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})]\right)$.

####  8.4.2 Zeller’s $g$-Prior

Similar to Section 8.3.1 (p. 8.3.1), suppose input $\bm{X}$ is transferred by
$\tilde{\bm{X}}=\bm{X}\bm{P}$ given some $p\times p$ matrix $\bm{P}$ in which
case, the model parameter is $\tilde{\boldsymbol{\beta}}$. Then, we have

$\bm{y}=\bm{X}\boldsymbol{\beta}=\tilde{\bm{X}}\tilde{\boldsymbol{\beta}}=\bm{X}\bm{P}\tilde{\boldsymbol{\beta}}.$

According to the principle of invariance, the posterior distributions of
$\boldsymbol{\beta}$ and $\bm{P}\tilde{\boldsymbol{\beta}}$ should be the
same. Simple calculation can show that this condition is met if
$\boldsymbol{\beta}_{0}=\boldsymbol{0},\boldsymbol{\Sigma}_{0}=g\sigma^{2}(\bm{X}^{\top}\bm{X})^{-1}$.
Following the Bayesian linear model with semi-conjugate prior, the posterior
of $\boldsymbol{\beta}$ is

$\mathrm{posterior}=\boldsymbol{\beta}|\bm{y},\bm{X},\sigma^{2}\sim\mathcal{N}(\boldsymbol{\beta}_{2},\boldsymbol{\Sigma}_{2}).$

where
$\boldsymbol{\Sigma}_{2}=\frac{g\sigma^{2}}{g+1}(\bm{X}^{\top}\bm{X})^{-1}$
and
$\boldsymbol{\beta}_{2}=\frac{g}{g+1}(\bm{X}^{\top}\bm{X})^{-1}(\bm{X}^{\top}\bm{y})$.

#### Derivation of $p(\bm{y}|\bm{X},\sigma^{2})$

Under the $g$-prior specified above, we derive the conditional distribution
$p(\bm{y}|\bm{X},\sigma^{2})$ which will be very useful for the Bayesian
variable selection procedure. We realize that

$\displaystyle p(\bm{y},\boldsymbol{\beta}|\bm{X},\sigma^{2})$
$\displaystyle=p(\bm{y}|\bm{X},\boldsymbol{\beta},\sigma^{2})p(\boldsymbol{\beta}|\bm{X},\sigma^{2})$
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{n/2}}\exp\left(-\frac{1}{2\sigma^{2}}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})\right)\qquad$
$\displaystyle(\text{$\bm{y}|\bm{X},\boldsymbol{\beta},\sigma^{2}\sim\mathcal{N}(\bm{X}\boldsymbol{\beta},\sigma^{2}\bm{I})$})$
$\displaystyle\,\,\,\,\,\,\times\frac{1}{(2\pi)^{n/2}|\boldsymbol{\Sigma}_{0}|^{1/2}}\exp\left(-\frac{1}{2}\boldsymbol{\beta}^{\top}\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}\right)\qquad$
$\displaystyle(\text{$\boldsymbol{\beta}|\bm{X},\sigma^{2}\sim\mathcal{N}(\boldsymbol{\beta}_{0},\boldsymbol{\Sigma}_{0})$})$
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{n/2}}\exp\left(-\frac{1}{2\sigma^{2}}\bm{y}^{\top}\bm{y}\right)\frac{|\boldsymbol{\Sigma}_{2}|^{1/2}}{|\boldsymbol{\Sigma}_{0}|^{1/2}}\exp\left(\frac{1}{2}\boldsymbol{\beta}_{2}^{\top}\boldsymbol{\Sigma}_{2}^{-1}\boldsymbol{\beta}_{2}\right)$
$\displaystyle\,\,\,\,\,\,\times\frac{1}{(2\pi)^{n/2}|\boldsymbol{\Sigma}_{2}|^{1/2}}\exp\left(-\frac{1}{2}(\boldsymbol{\beta}-\boldsymbol{\beta}_{2})^{\top}\boldsymbol{\Sigma}_{2}^{-1}(\boldsymbol{\beta}-\boldsymbol{\beta}_{2})\right),$

where $\boldsymbol{\beta}$ only appears in the third term. And the third term
is a multivariate normal with mean $\boldsymbol{\beta}_{2}$ and covariance
$\boldsymbol{\Sigma}_{2}$ and integrates to 1. And we notice that,
$\boldsymbol{\Sigma}_{2}=\frac{1}{g+1}\boldsymbol{\Sigma}_{0}$ such that
$\frac{|\boldsymbol{\Sigma}_{2}|^{1/2}}{|\boldsymbol{\Sigma}_{0}|^{1/2}}=\frac{1}{(g+1)^{p/2}}$.
Therefore, we obtain

$\displaystyle p(\bm{y}|\bm{X},\sigma^{2})$ $\displaystyle=\int
p(\bm{y},\boldsymbol{\beta}|\bm{X},\sigma^{2})d\boldsymbol{\beta}$ (8.4)
$\displaystyle=\int
p(\bm{y}|\bm{X},\boldsymbol{\beta},\sigma^{2})p(\boldsymbol{\beta}|\bm{X},\sigma^{2})d\boldsymbol{\beta}$
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{n/2}}\exp\left(-\frac{1}{2\sigma^{2}}\bm{y}^{\top}\bm{y}\right)\cdot\frac{1}{(g+1)^{p/2}}\exp\left(\frac{1}{2}\boldsymbol{\beta}_{2}^{\top}\boldsymbol{\Sigma}_{2}^{-1}\boldsymbol{\beta}_{2}\right)$
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{n/2}}\frac{1}{(g+1)^{p/2}}\exp\left(-\frac{r}{2\sigma^{2}}\right)$

where
$r=\bm{y}^{\top}\bm{y}-\bm{y}^{\top}(\frac{g}{g+1}\bm{X}(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top})\bm{y}$.

####  8.4.3 Bayesian Variable Selection

We previously introduce variable selection by $F$-test in Section 4.3 (p.
4.3). An alternative way can be achieved by Bayesian variable selection.

#### The Model

Suppose $\bm{z}=[z_{1},z_{2},\ldots,z_{p}]\in^{p}$ is a mask vector such that
$z_{j}\in\\{0,1\\}$ for all $j\in\\{1,2,\ldots,p\\}$. For each variable
$\beta_{j}$ in $\boldsymbol{\beta}$, we set $\beta_{j}=z_{j}\times b_{j}$
where $b_{j}$ can be understood as the original variable, and $\beta_{j}$ is
the final variable. That is $\boldsymbol{\beta}=\bm{z}\odot\bm{b}$ where
$\odot$ is the Hadamard product. Then the model can be expressed as

$\bm{y}=\bm{X}\boldsymbol{\beta}+\boldsymbol{\epsilon}=\bm{X}(\bm{z}\odot\bm{b})+\boldsymbol{\epsilon},$

where $\boldsymbol{\epsilon}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\bm{I})$.
Bayesian variable selection can be understood as obtaining a posterior
distribution for the mask vector $\bm{z}$. By Bayes’ theorem, we obtain the
posterior distribution by

$p(\bm{z}|\bm{y},\bm{X})\propto p(\bm{z})p(\bm{y}|\bm{X},\bm{z}).$

Alternatively, the ratio of two models $\bm{z}_{a}$, $\bm{z}_{b}$ is

$\displaystyle\text{odds}(\bm{z}_{a},\bm{z}_{b}|\bm{y},\bm{X})$
$\displaystyle=\frac{p(\bm{z}_{a}|\bm{y},\bm{X})}{p(\bm{z}_{b}|\bm{y},\bm{X})}$
$\displaystyle=$ $\displaystyle\,\,\,\,\frac{p(\bm{z}_{a})}{p(\bm{z}_{b})}$
$\displaystyle\times$
$\displaystyle\,\,\,\,\,\frac{p(\bm{y}|\bm{X},\bm{z}_{a})}{p(\bm{y}|\bm{X},\bm{z}_{b})}$
posterior odds $\displaystyle=$ prior odds $\displaystyle\times$ Bayes factor

where the Bayes factor can be understood as how much the data favor the model
$\bm{z}_{a}$ over the model $\bm{z}_{b}$.

#### Derivation of the Bayes Factor

Write out the Bayes factor by

$\displaystyle p(\bm{y}|\bm{X},\bm{z})$ $\displaystyle=\int\int
p(\bm{y},\boldsymbol{\beta},\sigma^{2}|\bm{X},\bm{z})d\boldsymbol{\beta}d\sigma^{2}$
$\displaystyle=\int\left(\int
p(\bm{y},\boldsymbol{\beta},|\bm{X},\bm{z},\sigma^{2})d\boldsymbol{\beta}\right)p(\sigma^{2})d\sigma^{2}$
$\displaystyle=\int
p(\bm{y}|\bm{X},\sigma^{2},\bm{z})p(\sigma^{2})d\sigma^{2},$

where $p(\bm{y}|\bm{X},\sigma^{2},\bm{z})=\left(\int
p(\bm{y},\boldsymbol{\beta},|\bm{X},\bm{z},\sigma^{2})d\boldsymbol{\beta}\right)$
can be obtained from Equation (8.4) (under Zeller’s $g$-prior) by substituting
$\bm{X}$ by $\bm{X}_{z}$ where we remove the variable $i$ if $z_{i}$=0.

We realize that $\gamma=\frac{1}{\sigma^{2}}$ and

$p(\sigma^{2})=p(\gamma|a_{0},b_{0})=\frac{{b_{0}}^{a_{0}}}{\Gamma(a_{0})}\gamma^{a_{0}-1}\exp(-b_{0}\gamma).$

Then,

$\displaystyle p(\bm{y}|\bm{X},\sigma^{2},\bm{z})p(\sigma^{2})$
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{n/2}}\frac{1}{(g+1)^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p_{z}}/2}}\exp\left(-\frac{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}r_{z}}}{2\sigma^{2}}\right)\cdot\frac{{b_{0}}^{a_{0}}}{\Gamma(a_{0})}\gamma^{a_{0}-1}\exp(-b_{0}\gamma)$
$\displaystyle=\frac{1}{(2\pi)^{n/2}}\frac{1}{(g+1)^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p_{z}}/2}}\frac{{b_{0}}^{a_{0}}}{\Gamma(a_{0})}\gamma^{a_{0}+n/2-1}\exp\left(-(b_{0}+\frac{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}r_{z}}}{2})\gamma\right)$
$\displaystyle=\frac{1}{(2\pi)^{n/2}}\frac{1}{(g+1)^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p_{z}}/2}}\frac{{b_{0}}^{a_{0}}}{\Gamma(a_{0})}\frac{\Gamma(a_{n})}{{b_{n}}^{a_{n}}}\cdot\frac{{b_{n}}^{a_{n}}}{\Gamma(a_{n})}\gamma^{a_{n}-1}\exp\left(-b_{n}\gamma\right)$
$\displaystyle=\frac{1}{(2\pi)^{n/2}}\frac{1}{(g+1)^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p_{z}}/2}}\frac{{b_{0}}^{a_{0}}}{\Gamma(a_{0})}\frac{\Gamma(a_{n})}{{b_{n}}^{a_{n}}}\cdot\mathrm{Gamma}(\gamma|a_{n},b_{n}),$

where
$r_{z}=\bm{y}^{\top}\bm{y}-\bm{y}^{\top}(\frac{g}{g+1}\bm{X}_{z}(\bm{X}_{z}^{\top}\bm{X}_{z})^{-1}\bm{X}_{z}^{\top})\bm{y}$,
$p_{z}$ is the number of 1’s in $\bm{z}$, $a_{n}=a_{0}+n/2$,
$b_{n}=b_{0}+\frac{r_{z}}{2}$, and $\mathrm{Gamma}(\gamma|a_{n},b_{n})$ is the
probability density function of Gamma distribution w.r.t. $\gamma$ with
parameters $a_{n},b_{n}$. Since $\gamma$ only appears in the last term
$\mathrm{Gamma}(\gamma|a_{n},b_{n})$ which integrates to 1, we have

$p(\bm{y}|\bm{X},\bm{z})=\frac{1}{(2\pi)^{n/2}}\frac{1}{(g+1)^{p_{z}/2}}\frac{{b_{0}}^{a_{0}}}{\Gamma(a_{0})}\frac{\Gamma(a_{n})}{{b_{n}}^{a_{n}}}.$

#### Same Prior Hyperparameters

Similarly, under models $\bm{z}_{a}$ and $\bm{z}_{b}$ where we assume the same
parameters of $a_{0}$ and $b_{0}$ in the two models, we have

$\frac{p(\bm{y}|\bm{X},\bm{z}_{a})}{p(\bm{y}|\bm{X},\bm{z}_{b})}=\frac{(g+1)^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p_{z_{b}}}/2}}{(g+1)^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p_{z_{a}}}/2}}\left(\frac{2b_{0}+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}r_{z_{b}}}}{2b_{0}+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}r_{z_{a}}}}\right)^{a_{0}+\frac{n}{2}},$

where $p_{z_{a}}$ is the number of variables selected in the model
$\bm{z}_{a}$, and $p_{z_{b}}$ is the number of variables selected in the model
$\bm{z}_{b}$.

#### Different Prior Hyperparameters

We have previously mentioned that

1\. $2a_{0}$ is the prior sample size for the noise
$\sigma^{2}=\frac{1}{\gamma}$.

2\. $\frac{2b_{0}}{a_{0}}$ is the prior variance of the noise.

Suppose now, given the two models $\bm{z}_{a}$ and $\bm{z}_{b}$, we assume
$2a_{0}=1$ for both of the models (i.e., prior sample sizes for the noise are
both 1), and set the $\frac{2b_{0}}{a_{0}}$ to be the estimated residual
variance under the least squares estimate for each model, say maximum
likelihood estimators
$\hat{\sigma}^{2}_{z_{a}}=\frac{1}{n}||\bm{y}-\bm{X}_{z_{a}}\hat{\boldsymbol{\beta}}_{z_{a}}||^{2}$
and
$\hat{\sigma}^{2}_{z_{b}}=\frac{1}{n}||\bm{y}-\bm{X}_{z_{b}}\hat{\boldsymbol{\beta}}_{z_{b}}||^{2}$
(which are biased estimators for $\sigma^{2}$, see Section 2.2, p. 2.2). Or we
could choose the unbiased estimators which are divided by $n-p_{z_{a}}$ and
$n-p_{z_{b}}$ respectively rather than divided by $n$ (see Section 3.4, p.
3.4). Then, we have

$\frac{p(\bm{y}|\bm{X},\bm{z}_{a})}{p(\bm{y}|\bm{X},\bm{z}_{b})}=\frac{(g+1)^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p_{z_{b}}}/2}}{(g+1)^{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}p_{z_{a}}}/2}}\left(\frac{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\hat{\sigma}^{2}_{z_{a}}}}{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\hat{\sigma}^{2}_{z_{b}}}}\right)^{\frac{1}{2}}\left(\frac{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\hat{\sigma}^{2}_{z_{b}}}+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}r_{z_{b}}}}{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\hat{\sigma}^{2}_{z_{a}}}+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}r_{z_{a}}}}\right)^{\frac{n+1}{2}}.$
(8.5)

Notice that the ratio of the marginal probabilities is essentially a balance
between the model complexity and goodness of fit:

  * •

A larger value of $p_{z_{b}}$ means the model $\bm{z}_{b}$ has more selected
variables (more complexity) which will make the ratio larger and penalize
model $\bm{z}_{b}$.

  * •

However, a more complex model will make $r_{z_{b}}$ smaller which in turn will
make the ratio smaller and penalize model $\bm{z}_{a}$.

#### Gibbs Sampler

Given a current value $\bm{z}=[z_{1},z_{2},\ldots,z_{p}]^{\top}$, a new value
of the $j$-th variable $z_{j}$ is generated by sampling from
$p(z_{j}|\bm{y},\bm{X},\bm{z}_{-j})$ where $z_{-j}$ refers to the values of
$\bm{z}$ except the $j$-th element $z_{j}$. Specifically, we define the
intermediate parameter222This intermediate parameter is quite useful in other
contexts, e.g., Bayesian inference for interpolative decomposition (Lu, 2022a,
b).

$o_{j}=\frac{p(z_{j}=1|\bm{y},\bm{X},\bm{z}_{-j})}{p(z_{j}=0|\bm{y},\bm{X},\bm{z}_{-j})}=\frac{p(z_{j}=1)}{p(z_{j}=0)}\times\frac{p(\bm{y}|\bm{X},\bm{z}_{-j},z_{j}=1)}{p(\bm{y}|\bm{X},\bm{z}_{-j},z_{j}=0)},$

where the last term can be obtained by Equation (8.5). Trivially, we can set
$p(z_{j}=1)=p(z_{j}=0)=0.5$. Then, using the intermediate parameter, the full
conditional probability of $z_{j}$ being equal to 1 can be obtained by

$p(z_{j}=1|\bm{y},\bm{X},\bm{z}_{-j})=\frac{o_{j}}{1+o_{j}}.$ (8.6)

Therefore, given the value of $\bm{z}^{(k)}$ in $k$-th step, we can generate
$\\{\bm{z}^{(k+1)},\gamma^{(k+1)},\boldsymbol{\beta}^{(k+1)}\\}$ by the
following steps:

0\. Set initial values to $\boldsymbol{\beta}$, $\gamma=\frac{1}{\sigma^{2}}$,
and $\bm{z}$ if $k$=1;

1\. Update $\bm{z}$: For $j\in\\{1,2,\ldots,p\\}$ in random order, replace
$z_{j}$ with a sample from $p(z_{j}=1|\bm{y},\bm{X},\bm{z}_{-j})$ (Equation
(8.6));

1\. Update $\boldsymbol{\beta}$:
$\boldsymbol{\beta}|\bm{y},\bm{X},\gamma,\bm{z}\sim\mathcal{N}(\boldsymbol{\beta}_{2},\boldsymbol{\Sigma}_{2})$,
where
$\boldsymbol{\Sigma}_{2}=\frac{g\sigma^{2}}{g+1}(\bm{X}_{z}^{\top}\bm{X}_{z})^{-1}$
and
$\boldsymbol{\beta}_{2}=\frac{g}{g+1}(\bm{X}_{z}^{\top}\bm{X}_{z})^{-1}(\bm{X}_{z}^{\top}\bm{y})$
(Equation (8.4.2), p. 8.4.2);

2\. Update $\gamma$:
$\gamma|\bm{y},\bm{X},\boldsymbol{\beta},\bm{z}\sim\mathrm{Gamma}\left(a_{0}+\frac{n}{2},[b_{0}+\frac{1}{2}||\bm{y}-\bm{X}_{z}\boldsymbol{\beta}_{z}||^{2}]\right)$
(Equation (8.1), p. 8.1).

###  8.5 Bayesian Linear Model: Full Conjugate Prior

To derive the full conjugate algorithm for the linear model, we place an
inverse-gamma prior over the variance parameter instead. Putting a gamma prior
over the inverse variance is equivalent to putting an inverse-gamma prior on
the variance.

######  Definition 8.4 (Inverse-Gamma Distribution)

A random variable $Y$ is said to follow the inverse-gamma distribution with
parameters $r>0$ and $\lambda>0$ if

$f(y;r,\lambda)=\left\\{\begin{aligned}
&\frac{\lambda^{r}}{\Gamma(r)}y^{-r-1}\exp(-\frac{\lambda}{y}),&\mathrm{\,\,if\,\,}y>0.\\\
&0,&\mathrm{\,\,if\,\,}y\leq 0,\end{aligned}\right.$

where again $\Gamma(\cdot)$ is the gamma function (Definition 3.1, p. 3.1).
And it is denoted by $Y\sim\mathrm{Inverse\mathchar 45\relax
Gamma}(r,\lambda)$. The mean and variance of inverse-gamma distribution are
given by

$\mathrm{E}[Y]=\left\\{\begin{aligned}
&\frac{\lambda}{r-1},\,&\mathrm{if\,}r\geq 1.\\\
&\infty,\,&\mathrm{if\,}0<r<1.\end{aligned}\right.\qquad\mathrm{Var}[Y]=\left\\{\begin{aligned}
&\frac{\lambda^{2}}{(r-1)^{2}(r-2)},\,&\mathrm{if\,}r\geq 2.\\\
&\infty,\,&\mathrm{if\,}0<r<2.\end{aligned}\right.$

Figure 8.4 compares different parameters for gamma distribution and inverse-
gamma distribution.

(a) Inverse-amma probability density functions for different values of the
parameters $r,\lambda$.

(b) Gamma probability density functions for different values of the parameters
$r,\lambda$.

Figure 8.4: Comparison between the inverse-gamma distribution and gamma
distribution (Definition 8.3, p. 8.3).

Note that the inverse-gamma density is not simply the gamma density with $x$
replaced by $\frac{1}{y}$. There is an additional factor of $y^{-2}$ from the
Jacobian in the change-of-variables formula. 333Which is from the Jacobian in
the change-of-variables formula. A short proof is provided here. Let
$y=\frac{1}{x}$ where $y\sim\mathrm{Inverse\mathchar 45\relax
Gamma}(r,\lambda)$ and $x\sim\mathrm{Gamma}(r,\lambda)$. Then,
$f(y)|dy|=f(x)|dx|$ which results in
$f(y)=f(x)|\frac{dx}{dy}|=f(x)x^{2}\xlongequal{\mathrm{y}=\frac{1}{x}}\frac{\lambda^{r}}{\Gamma(r)}y^{-r-1}exp(-\frac{\lambda}{y})$
for $y>0$.

Same setting as the semiconjugate prior distribution in Section 8.4 (p. 8.4).
We have the likelihood function:

$\mathrm{likelihood}=\bm{y}|\bm{X},\boldsymbol{\beta},\sigma^{2}\sim\mathcal{N}(\bm{X}\boldsymbol{\beta},\sigma^{2}\bm{I}),$

which is the same as the likelihood density in the zero-mean model (Section
8.3, p. 8.3) and the semiconjugate model (Section 8.4, p. 8.4). But now we
specify a Gaussian prior (with unfixed covariance matrix) over the weight
parameter, and an inverse-gamma density over the variance parameter by

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathrm{prior:\,}}$
$\displaystyle\boldsymbol{\beta}|\sigma^{2}\sim\mathcal{N}(\boldsymbol{\beta}_{0},{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\sigma^{2}}\boldsymbol{\Sigma}_{0})$
$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\sigma^{2}\sim\mathrm{Inverse\mathchar
45\relax Gamma}(a_{0},b_{0})},$

where again we differentiate from previous descriptions by blue text. Note we
place an inverse-gamma density over the variance parameter in the full
conjugate model instead of a gamma density over the precision (inverse
variance) parameter in the semiconjugate case. However, it is demonstrable
that the two circumstances are identical (via the Jacobian in the change-of-
variables formula).

Equivalently, we can formulate the prior into one which is called the normal-
inverse-gamma (NIG) distribution:

$\displaystyle\mathrm{prior:\,}$
$\displaystyle\boldsymbol{\beta},\sigma^{2}\sim
NIG(\boldsymbol{\beta}_{0},\boldsymbol{\Sigma}_{0},a_{0},b_{0})$
$\displaystyle=\mathcal{N}(\boldsymbol{\beta}_{0},\sigma^{2}\boldsymbol{\Sigma}_{0})\cdot\mathrm{Inverse\mathchar
45\relax Gamma}(a_{0},b_{0}).$

Again by applying the Bayes’ theorem
“$\mathrm{posterior}\propto\mathrm{likelihood}\times\mathrm{prior}$”, we
obtain the posterior

$\displaystyle\mathrm{posterior}$
$\displaystyle=p(\boldsymbol{\beta},\sigma^{2}|\bm{y},\bm{X})$
$\displaystyle\propto
p(\bm{y}|\bm{X},\boldsymbol{\beta},\sigma^{2})p(\boldsymbol{\beta},\sigma^{2}|\boldsymbol{\beta}_{0},\boldsymbol{\Sigma}_{0},a_{0},b_{0})$
$\displaystyle=\frac{1}{(2\pi\sigma^{2})^{n/2}}\exp\left\\{-\frac{1}{2\sigma^{2}}(\bm{y}-\bm{X}\boldsymbol{\beta})^{\top}(\bm{y}-\bm{X}\boldsymbol{\beta})\right\\}$
$\displaystyle\,\,\,\,\,\,\times\frac{1}{(2\pi\sigma^{2})^{p/2}|\boldsymbol{\Sigma}_{0}|^{1/2}}\exp\left\\{-\frac{1}{2\sigma^{2}}(\boldsymbol{\beta}-\boldsymbol{\beta}_{0})^{\top}\boldsymbol{\Sigma}_{0}^{-1}(\boldsymbol{\beta}-\boldsymbol{\beta}_{0})\right\\}$
$\displaystyle\,\,\,\,\,\,\times\frac{{b_{0}}^{a_{0}}}{\Gamma(a_{0})}\frac{1}{(\sigma^{2})^{a_{0}+1}}\exp\\{-\frac{b_{0}}{\sigma^{2}}\\}$
$\displaystyle\propto\frac{1}{(2\pi\sigma^{2})^{p/2}}\exp\left\\{\frac{1}{2\sigma^{2}}(\boldsymbol{\beta}-\boldsymbol{\beta}_{3})^{\top}\boldsymbol{\Sigma}_{3}^{-1}(\boldsymbol{\beta}-\boldsymbol{\beta}_{3})\right\\}$
$\displaystyle\,\,\,\,\,\,\times\frac{1}{(\sigma^{2})^{a_{0}+\frac{n}{2}+1}}\exp\left\\{-\frac{1}{\sigma^{2}}[b_{0}+\frac{1}{2}(\bm{y}^{\top}\bm{y}+\boldsymbol{\beta}_{0}^{\top}\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0}-\boldsymbol{\beta}_{3}^{\top}\boldsymbol{\Sigma}_{3}^{-1}\boldsymbol{\beta}_{3})]\right\\},$

where
$\boldsymbol{\Sigma}_{3}=(\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}$
and

$\boldsymbol{\beta}_{3}=\boldsymbol{\Sigma}_{3}(\bm{X}^{\top}\bm{y}+\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0})=(\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}(\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0}+\bm{X}^{\top}\bm{y}).$

Let $a_{n}=a_{0}+\frac{n}{2}+1$ and
$b_{n}=b_{0}+\frac{1}{2}(\bm{y}^{\top}\bm{y}+\boldsymbol{\beta}_{0}^{\top}\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0}-\boldsymbol{\beta}_{3}^{\top}\boldsymbol{\Sigma}_{3}^{-1}\boldsymbol{\beta}_{3})$.
The posterior admits conjugacy and is a NIG distribution:

$\displaystyle\mathrm{posterior}$
$\displaystyle=\boldsymbol{\beta},\sigma^{2}|\bm{y},\bm{X}\sim
NIG(\boldsymbol{\beta}_{3},\boldsymbol{\Sigma}_{3},a_{n},b_{n}).$

1\. $\boldsymbol{\Sigma}_{0}$ here is a fixed hyperparameter. 2\. If we assume
further $\bm{X}$ has full rank, when
$\boldsymbol{\Sigma}_{0}^{-1}\rightarrow\boldsymbol{0}$,
$\boldsymbol{\beta}_{3}\rightarrow\hat{\boldsymbol{\beta}}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}\bm{y}$
which is the OLS estimator. 3\. When $b_{0}\rightarrow\infty$, then
$\sigma^{2}\rightarrow\infty$ and $\boldsymbol{\beta}_{3}$ is approximately
$\boldsymbol{\beta}_{0}$, the prior expectation of parameter. Compared to
$\boldsymbol{\beta}_{2}$ in Section 8.4 (p. 8.4),
$\sigma^{2}\rightarrow\infty$ will make $\boldsymbol{\beta}_{2}$ approach to
$\boldsymbol{\beta}_{0}$ where $\sigma^{2}$ is a fixed hyperparameter. 4\.
Weighted average: we reformulate $\displaystyle\boldsymbol{\beta}_{3}$
$\displaystyle=(\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}(\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0}+\bm{X}^{\top}\bm{y})$
$\displaystyle=(\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}\boldsymbol{\Sigma}_{0}^{-1}\boldsymbol{\beta}_{0}+(\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}(\bm{X}^{\top}\bm{X})(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{y}$
$\displaystyle=(\bm{I}-\bm{C})\boldsymbol{\beta}_{0}+\bm{C}\hat{\boldsymbol{\beta}},$
where $\hat{\boldsymbol{\beta}}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{y}$
is the OLS estimator of $\boldsymbol{\beta}$ and
$\bm{C}=(\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}(\bm{X}^{\top}\bm{X})$.
We see that the posterior mean of $\boldsymbol{\beta}$ is a weighted average
of the prior mean and the OLS estimator of $\boldsymbol{\beta}$. Thus, if we
set $\boldsymbol{\beta}_{0}=\hat{\boldsymbol{\beta}}$, the posterior mean of
$\boldsymbol{\beta}$ will be exactly $\hat{\boldsymbol{\beta}}$. 5\. From
$a_{n}=a_{0}+\frac{n}{2}+1$, we realize that $2a_{0}$ is the prior sample size
for $\sigma^{2}$. 6\.
$\boldsymbol{\Sigma}_{3}^{-1}=\bm{X}^{\top}\bm{X}+\boldsymbol{\Sigma}_{0}^{-1}$:
The posterior inverse covariance is equal to $\bm{X}^{\top}\bm{X}$ \+ prior
inverse covariance.

##  Chapter 9 Beyond Bayesian Approach: Gaussian Process Regression

###  9.1 Beyond Bayesian Approach: Gaussian Process Regression

In some textbooks or college courses, people talk about the Gaussian process
without any introduction to the Bayesian linear model. This is quite
confusing. We will discover that the Bayesian linear model with zero-mean
prior is where the Gaussian process regression genuinely ”arises.”

###  9.2 Predictive Distribution of Bayesian Linear Model with Zero-Mean Prior

Following from Bayesian linear model with zero mean in Section 8.3 (p. 8.3).
The predictive distribution $g_{\ast}=g_{\ast}(\bm{x}_{\ast})$ for a new data
point $\bm{x}_{\ast}$ is again a Gaussian distribution

$\displaystyle p(g_{\ast}|\bm{x}_{\ast}\bm{X},\bm{y},\sigma^{2})$
$\displaystyle=\int
p(g_{\ast}|\bm{x}_{\ast},\boldsymbol{\beta})p(\boldsymbol{\beta}|\bm{X},\bm{y},\sigma^{2})d\boldsymbol{\beta}$
(9.1)
$\displaystyle=\mathcal{N}(\frac{1}{\sigma^{2}}\bm{x}_{\ast}^{\top}\boldsymbol{\Sigma}_{1}\bm{X}^{\top}\bm{y},\bm{x}_{\ast}^{\top}\boldsymbol{\Sigma}_{1}\bm{x}_{\ast})$
$\displaystyle=\mathcal{N}(\bm{x}_{\ast}^{\top}\boldsymbol{\beta}_{1},\bm{x}_{\ast}^{\top}\boldsymbol{\Sigma}_{1}\bm{x}_{\ast}),$

where the mean of the predictive distribution is the posterior mean of weight
parameter (i.e., $\boldsymbol{\beta}_{1}$) multiplied by the new input. And
the predictive variance is a quadratic form of the new input indicating that
the predictive uncertainties grow with the magnitude of the new input.

However, this Bayesian linear model suffers from limited expressiveness. One
solution to overcome this issue is to use a set of basis functions and
transfer from the $p\mathchar 45\relax$dimensional space to higher dimensional
space, say $q\mathchar 45\relax$dimensional space:
$\bm{x}\in^{p}\rightarrow\boldsymbol{\phi}(\bm{x})\in^{q}$ and
$\bm{X}\rightarrow\boldsymbol{\Phi}(\bm{X})\in^{n\times q}$. But another issue
arises that the computation grows quadratically, e.g.,
$\mathcal{O}(p^{2}+p)\rightarrow\mathcal{O}(q^{2}+q)$ for the computation of
predictive variance.

######  Remark 9.1 (Kernel Trick)

Kernel trick could help reduce computational complexity in high dimensional
space. The key point to use kernel trick is to make sure that the input space
only involves dot products.

######  Example 9.2.1 (Kernel Trick)

1\. The computation of
$\boldsymbol{\phi}(\bm{x}_{\ast})^{\top}\boldsymbol{\Phi}(\bm{X})^{\top}=\bm{z}\in^{1\times
n}$ only has the dot product from the input space. So we could use kernel
trick to do the computation. And each $z_{i}=k(\bm{x}_{\ast},\bm{x}_{i})$
where $k(\cdot,\cdot)$ is a kernel function. The computation of kernel
function is potentially less than $\mathcal{O}(q)$. 2\. The computation of
$\bm{y}^{\top}\boldsymbol{\Phi}(\bm{X})$ has dot product between the input
space and output space. So it cannot be reduced to kernel form.

Having this in mind, if we can reduce the computation of Equation (9.1)
involving input space to having only dot products, we could apply the kernel
trick.

Transfer the parameters in Equation (9.1) into high dimensional space. By
transforming from {$\bm{x},\bm{X}$} in $x$-space $\in^{p}$ into
{$\boldsymbol{\phi},\boldsymbol{\Phi}$} in $z$-space $\in^{q}$, we obtain the
predictive distribution in the $z$-space form:

$\boxed{\begin{aligned} z\mathchar 45\relax\mathrm{Space\,Form:}\\\
g_{\ast}|\bm{x}_{\ast}\bm{X},\bm{y},\sigma^{2}&\sim\mathcal{N}(\boldsymbol{\phi}(\bm{x}_{\ast})^{\top}\boldsymbol{\beta}_{1},\boldsymbol{\phi}(\bm{x}_{\ast})^{\top}\boldsymbol{\Sigma}_{1}\boldsymbol{\phi}(\bm{x}_{\ast}))\\\
&=\mathcal{N}\left(\boldsymbol{\phi}(\bm{x}_{\ast})^{\top}(\frac{1}{\sigma^{2}}\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}(\frac{1}{\sigma^{2}}\boldsymbol{\Phi}^{\top}\bm{y}),\boldsymbol{\phi}(\bm{x}_{\ast})^{\top}\boldsymbol{\Sigma}_{1}\boldsymbol{\phi}(\bm{x}_{\ast})\right),\end{aligned}}$

where
$\boldsymbol{\Sigma}_{1}=(\frac{1}{\sigma^{2}}\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\boldsymbol{\Sigma}_{0}^{-1})^{-1}$.
Let $\bm{K}=\boldsymbol{\Phi}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top}$
and $\boldsymbol{\phi}_{\ast}=\boldsymbol{\phi}(\bm{x}_{\ast})$, we have

$\frac{1}{\sigma^{2}}\boldsymbol{\Phi}^{\top}(\bm{K}+\sigma^{2}\bm{I})=\frac{1}{\sigma^{2}}\boldsymbol{\Phi}^{\top}(\boldsymbol{\Phi}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top}+\sigma^{2}\bm{I})=\boldsymbol{\Sigma}_{1}^{-1}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top},$
(9.2)

where we use the fact that
$\bm{A}\bm{B}\bm{A}+\bm{A}=\bm{A}(\bm{B}\bm{A}+\bm{I})=(\bm{A}\bm{B}+\bm{I})\bm{A}$.
Note that we do not differentiate the notations of $\boldsymbol{\Sigma}_{0}$,
$\boldsymbol{\Sigma}_{1}$ in $x$ and $z$-spaces. But only differentiate the
notation of $\bm{X}$ and $\boldsymbol{\Phi}$ in $x$ and $z$-spaces
respectively. Then by multiplying Equation (9.2) by $\boldsymbol{\Sigma}_{1}$
from left and $(\bm{K}+\sigma^{2}\bm{I})^{-1}$ from right, we transform the
predictive prediction into inner product form:

$\boxed{\begin{aligned} \mathrm{Inner\,Product\,Form:}\\\
g_{\ast}|\bm{x}_{\ast}\bm{X},\bm{y},\sigma^{2}\sim\mathcal{N}(&\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top}(\bm{K}+\sigma^{2}\bm{I})^{-1}\bm{y},\\\
&\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}_{\ast}-\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top}(\bm{K}+\sigma^{2}\bm{I})^{-1}\boldsymbol{\Phi}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}_{\ast}),\end{aligned}}$

where the input space only has the forms
$\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}_{\ast},\,\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top},\,\boldsymbol{\Phi}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}$.
So we meet the kernel trick requirement.

Define the kernel variables:
$K(\bm{x}_{\ast},\bm{x}_{\ast})=\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}_{\ast}$,
$K(\bm{x}_{\ast},\bm{X})=\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top}$,
$K(\bm{X},\bm{X})=\boldsymbol{\Phi}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top}$.
Then we transform the predictive distribution into the kernel form:

$\boxed{\begin{aligned} \mathrm{Kernel\,Form:}\\\
g_{\ast}|\bm{x}_{\ast}\bm{X},\bm{y},\sigma^{2}\sim\mathcal{N}(&K(\bm{x}_{\ast},\bm{X})(K(\bm{X},\bm{X})+\sigma^{2}\bm{I})^{-1}\bm{y},\\\
&K(\bm{x}_{\ast},\bm{x}_{\ast})-K(\bm{x}_{\ast},\bm{X})(K(\bm{X},\bm{X})+\sigma^{2}\bm{I})^{-1}K(\bm{X},\bm{x}_{\ast})).\end{aligned}}$
(9.3)

###  9.3 Kernels In A Nutshell

Kernel with basis functions transforms
$\bm{x}\in^{p}\rightarrow\boldsymbol{\phi}(\bm{x})\in^{q}$. Thus the inner
product in the $p\mathchar 45\relax$dimensional space
$\bm{x}^{\top}\bm{x}^{\prime}$ is transformed into the inner product in the
$q\mathchar 45\relax$dimensional space
$k(\bm{x},\bm{x}^{\prime})=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\phi}(\bm{x}^{\prime})$.
Then the kernel matrix $K(\bm{X},\bm{X})$ has the following properties

1\. $K$ should be symmetric, i.e.,
$k(\bm{x},\bm{x}^{\prime})=k(\bm{x}^{\prime},\bm{x})$.

2\. The kernel matrix $K$ is positive semi-definite (PSD).

Proof [of Kernel Matrix $K$ is PSD] Let $K_{ij}=k(\bm{x}_{i},\bm{x}_{j})$,
$\forall i,j\in\\{1,2,...,n\\}$. And for any vector $\bm{t}$, we have

$\displaystyle\bm{t}^{\top}K\bm{t}$
$\displaystyle=\sum_{i,j=1}^{n}t_{i}t_{j}K_{ij}$
$\displaystyle=\sum_{i,j=1}^{n}t_{i}t_{j}\boldsymbol{\phi}(\bm{x}_{i})^{\top}\boldsymbol{\phi}(\bm{x}_{j})$
$\displaystyle=(\sum_{i=1}^{n}t_{i}\boldsymbol{\phi}(\bm{x}_{i}))^{\top}(\sum_{j=1}^{n}t_{j}\boldsymbol{\phi}(\bm{x}_{j}))$
$\displaystyle=||\sum_{i=1}^{n}t_{i}\boldsymbol{\phi}(\bm{x}_{i})||^{2}\geq
0.$

This completes the proof.

We may find $k(\bm{x},\bm{x}^{\prime})$ can be any function related to
$\bm{x},\bm{x}^{\prime}$ at a first glance. But from this PSD property of the
kernel matrix, this function is restricted to a small part in order to make
the $k(\bm{x},\bm{x}^{\prime})$ to be an inner-product form implicitly in
$q$-dimensional space.

######  Remark 9.2 (Some Specific Kernels)

1\. Linear Kernel: $k(\bm{x},\bm{x}^{\prime})=\bm{x}^{\top}\bm{x}^{\prime}$.
2\. Polynomial Kernel:
$k(\bm{x},\bm{x}^{\prime})=(\eta+\gamma\bm{x}^{\top}\bm{x}^{\prime})^{Q}$ with
$\gamma>0,\eta\geq 0$. 3\. Gaussian Kernel:
$k(\bm{x},\bm{x}^{\prime})=\exp(-\gamma||\bm{x}-\bm{x}^{\prime}||^{2})$. We
here prove that Gaussian kernel is of infinite-dimensional transformation.
Without loss of generality, let $\gamma=1$, then, $\displaystyle
k(\bm{x},\bm{x}^{\prime})$
$\displaystyle=\exp(-||\bm{x}-\bm{x}^{\prime}||^{2})$
$\displaystyle=\exp(-\bm{x}^{\top}\bm{x})\exp(-\bm{x}^{\prime\top}\bm{x}^{\prime})\exp(2\bm{x}^{\top}\bm{x}^{\prime})$
$\displaystyle\underset{\mathrm{expansion}}{\overset{\mathrm{Taylor}}{=}}\exp(-\bm{x}^{\top}\bm{x})\exp(-\bm{x}^{\prime\top}\bm{x}^{\prime})\exp\left(\sum_{i=0}^{\infty}\frac{(2\bm{x}^{\top}\bm{x}^{\prime})^{i}}{i!}\right)$
$\displaystyle=\sum_{i=0}^{\infty}\left(\exp(-\bm{x}^{\top}\bm{x})\exp(-\bm{x}^{\prime\top}\bm{x}^{\prime})\sqrt{\frac{2^{i}}{i!}}\sqrt{\frac{2^{i}}{i!}}(\bm{x})^{i}\cdot(\bm{x}^{\prime})^{i}\right)$
$\displaystyle=\sum_{i=0}^{\infty}\left({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\exp(-\bm{x}^{\top}\bm{x})\sqrt{\frac{2^{i}}{i!}}(\bm{x})^{i}}\cdot{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\exp(-\bm{x}^{\prime\top}\bm{x}^{\prime})\sqrt{\frac{2^{i}}{i!}}(\bm{x}^{\prime})^{i}}\right)$
$\displaystyle={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\boldsymbol{\phi}(\bm{x})^{\top}}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{\phi}(\bm{x}^{\prime})},$
where
$\boldsymbol{\phi}(\bm{x})=\sum_{i=0}^{\infty}\exp(-\bm{x}^{\top}\bm{x})\sqrt{\frac{2^{i}}{i!}}(\bm{x})^{i}$,
which is a sum of 1-dimensional input to infinite-dimensional input.
Similarly, we can get the infinite-dimensional transformation of Gaussian
kernel when $\gamma\neq 1$. 4\. Other Valid Kernels: A big advantage of using
kernels is that we do not need to specify $\boldsymbol{\phi}(\bm{x})$
explicitly, since we can work directly with $K$.

###  9.4 Gaussian Process from Zero-Mean Prior

We use Gaussian processes (GPs) to describe a distribution over functions.
Gaussian processes are natural generalizations of multivariate Gaussian random
variables to infinite (countably or continuous) index sets. Formally we define
Gaussian processes as follows:

######  Definition 9.1 (Gaussian Process)

A Gaussian process is a collection of random variables, any finite number of
which have a joint Gaussian distribution. The definition does not exclude
Gaussian processes with finite index sets which would be simply Gaussian
distributions.

####  9.4.1 Noise-Free Observations

Following the Bayesian linear model with zero-mean prior in Section 8.3 (p.
8.3), we assume zero-mean prior
$\boldsymbol{\beta}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{0})$
and define $g(\bm{x})=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\beta}$ for
each input observation $\bm{x}$. Then the mean and covariance for the prior
output are:

$\displaystyle\mathrm{E}[g(\bm{x})]$
$\displaystyle=\boldsymbol{\phi}(\bm{x})^{\top}\mathrm{E}[\boldsymbol{\beta}]=0,$
(9.4) $\displaystyle\mathrm{E}[g(\bm{x})g(\bm{x}^{\prime})]$
$\displaystyle=\boldsymbol{\phi}(\bm{x})^{\top}\mathrm{E}[\boldsymbol{\beta}\boldsymbol{\beta}^{\top}]\boldsymbol{\phi}(\bm{x}^{\prime})=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}(\bm{x}^{\prime}),$

where the prior covariance matrix is specified manually and we could use a
kernel function instead:
$k(\bm{x},\bm{x}^{\prime})=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}(\bm{x}^{\prime})$.

Suppose we have the training input design matrix $\bm{X}$, training output
vector $\bm{y}$, test input design matrix $\bm{X}_{\ast}$, and test output
vector $\bm{g}_{\ast}$. We obtain the joint distribution of the training
outputs $\bm{y}$ and the test outputs $\boldsymbol{g}_{\ast}$ by Equation
(9.4):

$\left[\begin{matrix}\bm{y}\\\
\boldsymbol{g}_{\ast}\end{matrix}\right]\sim\mathcal{N}\left(\boldsymbol{0},\left[\begin{matrix}K(\bm{X},\bm{X}),&K(\bm{X},\boldsymbol{X}_{\ast})\\\
K(\boldsymbol{X}_{\ast},\bm{X}),&K(\boldsymbol{X}_{\ast},\boldsymbol{X}_{\ast})\end{matrix}\right]\right).$

Therefore, we could use Gaussian identities to get the marginal distribution
of test outputs $\boldsymbol{g}_{\ast}$.

######  Lemma 9.3 (Marginal Distribution of Test Outputs)

Follow from Section 8.3 (p. 8.3), assume zero-mean prior
$\boldsymbol{\beta}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{0})$
and define $g(\bm{x})=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\beta}$.
Given observed training inputs $\bm{X}$, training outputs $\bm{y}$ and test
inputs $\bm{x}_{\ast}$, the marginal distribution of test outputs
$\boldsymbol{g}_{\ast}$ is

$\displaystyle\bm{g}_{\ast}|\bm{X}_{\ast}\bm{X},\bm{y}\sim\mathcal{N}($
$\displaystyle K(\bm{X}_{\ast},\bm{X})K(\bm{X},\bm{X})^{-1}\bm{y},$
$\displaystyle
K(\bm{X}_{\ast},\bm{x}_{\ast})-K(\bm{X}_{\ast},\bm{X})K(\bm{X},\bm{X})^{-1}K(\bm{X},\bm{X}_{\ast})).$

Proof [of Lemma 9.3] Let $\bm{x}$ and $\bm{y}$ be jointly Gaussian random
vectors:

$\left[\begin{matrix}\bm{x}\\\
\bm{y}\end{matrix}\right]\sim\mathcal{N}\left(\left[\begin{matrix}\bm{u}_{x}\\\
\bm{u}_{y}\end{matrix}\right],\left[\begin{matrix}\bm{A},&\bm{C}\\\
\bm{C}^{\top},&\bm{B}\end{matrix}\right]\right).$ By Gaussian identity, the
marginal distribution of $\bm{x}$ and the conditional distribution of $\bm{x}$
given $\bm{y}$ are $\displaystyle\bm{x}$
$\displaystyle\sim\mathcal{N}(\bm{u}_{x},\bm{A}),$
$\displaystyle\bm{x}|\bm{y}$
$\displaystyle\sim\mathcal{N}(\bm{u}_{x}+\bm{C}\bm{B}^{-1}(\bm{y}-\bm{u}_{y}),\bm{A}-\bm{C}\bm{B}^{-1}\bm{C}^{\top})$
$\displaystyle\bm{y}$ $\displaystyle\sim\mathcal{N}(\bm{u}_{y},\bm{B})$
$\displaystyle\bm{y}|\bm{x}$
$\displaystyle\sim\mathcal{N}(\bm{u}_{y}+\bm{C}^{\top}\bm{A}^{-1}(\bm{x}-\bm{u}_{x}),\bm{B}-\bm{C}^{\top}\bm{A}^{-1}\bm{C}).$

Thus,

$\displaystyle\bm{g}_{\ast}|\bm{X}_{\ast}\bm{X},\bm{y}\sim\mathcal{N}($
$\displaystyle K(\bm{X}_{\ast},\bm{X})K(\bm{X},\bm{X})^{-1}\bm{y},$
$\displaystyle
K(\bm{X}_{\ast},\bm{x}_{\ast})-K(\bm{X}_{\ast},\bm{X})K(\bm{X},\bm{X})^{-1}K(\bm{X},\bm{X}_{\ast})),$

which completes the proof.

The marginal distribution of test outputs $\boldsymbol{g}_{\ast}$ agrees with
the kernel form of the Bayesian linear model with zero-mean prior in Equation
(9.3) except the noise term and the output are vectors now (i.e.,
$g_{\ast}\rightarrow\bm{g}_{\ast}$).

When given
$k(\bm{x},\bm{x}^{\prime})=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}(\bm{x}^{\prime})$
which implies a zero-mean prior
$\boldsymbol{\beta}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{0})$,
we have prior distribution over $\boldsymbol{g}_{\ast}$
$\begin{matrix}\boldsymbol{g}_{\ast}\end{matrix}\sim\mathcal{N}\left(\boldsymbol{0},\left[\begin{matrix}K(\boldsymbol{X}_{\ast},\boldsymbol{X}_{\ast})\end{matrix}\right]\right).$
Thus we can sample functions from
$p(\boldsymbol{g}_{\ast}|\boldsymbol{X}_{\ast})$. By sampling functions, we
mean the realizations of the output values given the inputs (possibly finite
or infinite number of samples) and the distribution of it. And this is the
meaning of distribution over functions. When given observed data $\bm{y}$ and
$\bm{X}$, we have prior distribution over $\bm{y},\boldsymbol{g}_{\ast}$
$\left[\begin{matrix}\bm{y}\\\
\boldsymbol{g}_{\ast}\end{matrix}\right]\sim\mathcal{N}\left(\boldsymbol{0},\left[\begin{matrix}K(\bm{X},\bm{X}),&K(\bm{X},\boldsymbol{X}_{\ast})\\\
K(\boldsymbol{X}_{\ast},\bm{X}),&K(\boldsymbol{X}_{\ast},\boldsymbol{X}_{\ast})\end{matrix}\right]\right).$
Thus we can compute the posterior by
$\displaystyle\bm{g}_{\ast}|\bm{X}_{\ast}\bm{X},\bm{y}\sim\mathcal{N}($
$\displaystyle K(\bm{X}_{\ast},\bm{X})K(\bm{X},\bm{X})^{-1}\bm{y},$
$\displaystyle
K(\bm{X}_{\ast},\bm{x}_{\ast})-K(\bm{X}_{\ast},\bm{X})K(\bm{X},\bm{X})^{-1}K(\bm{X},\bm{X}_{\ast})),$
by which we can sample functions for the posterior distribution of
$p(\bm{g}_{\ast}|\bm{X}_{\ast}\bm{X},\bm{y})$.

####  9.4.2 Noisy Observations

Following again from Section 8.3 (p. 8.3]), we assume zero-mean prior
$\boldsymbol{\beta}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{0})$.
And we define
$g(\bm{x})=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\beta}+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{\epsilon}}$
and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{\epsilon}\sim\mathcal{N}(0,\sigma^{2})}$.
Then the mean and covariance for the prior output are:

$\displaystyle\mathrm{E}[g(\bm{x})]$
$\displaystyle=\boldsymbol{\phi}(\bm{x})^{\top}\mathrm{E}[\boldsymbol{\beta}]+\mathrm{E}[\boldsymbol{\epsilon}]=\boldsymbol{0},$
$\displaystyle\mathrm{E}[g(\bm{x})g(\bm{x})]$
$\displaystyle=\boldsymbol{\phi}(\bm{x})^{\top}\mathrm{E}[\boldsymbol{\beta}\boldsymbol{\beta}^{\top}]\boldsymbol{\phi}(\bm{x})+2\boldsymbol{\phi}(\bm{x})^{\top}\mathrm{E}[\boldsymbol{\beta}]\mathrm{E}[\epsilon]+\mathrm{E}[\epsilon^{2}]$
$\displaystyle=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}(\bm{x})+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\sigma^{2}},$
$\displaystyle\mathrm{E}[g(\bm{x})g(\bm{x}^{\prime})]$
$\displaystyle=\boldsymbol{\phi}(\bm{x})^{\top}\mathrm{E}[\boldsymbol{\beta}\boldsymbol{\beta}^{\top}]\boldsymbol{\phi}(\bm{x}^{\prime})+\boldsymbol{\phi}(\bm{x})^{\top}\mathrm{E}[\boldsymbol{\beta}]\mathrm{E}[\epsilon]+\boldsymbol{\phi}(\bm{x}^{\prime})^{\top}\mathrm{E}[\boldsymbol{\beta}]\mathrm{E}[\epsilon]++\mathrm{E}[\epsilon^{2}]$
$\displaystyle=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}(\bm{x}^{\prime}).$

That is

$\mathrm{cov}(y_{i},y_{j})=k(\bm{x}_{i},\bm{x}_{j})\delta_{ij},$

where $\delta_{ij}$ is a Kronecker delta which is equal to 1 if and only if
when $i=j$ and 0 otherwise. Thus, we obtain the joint distribution of the
training outputs $\bm{y}$ and the test outputs $\boldsymbol{g}_{\ast}$

$\left[\begin{matrix}\bm{y}\\\
\boldsymbol{g}_{\ast}\end{matrix}\right]\sim\mathcal{N}\left(\boldsymbol{0},\left[\begin{matrix}K(\bm{X},\bm{X})+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\sigma^{2}\bm{I}},&K(\bm{X},\boldsymbol{X}_{\ast})\\\
K(\boldsymbol{X}_{\ast},\bm{X}),&K(\boldsymbol{X}_{\ast},\boldsymbol{X}_{\ast})\end{matrix}\right]\right).$

Similarly, we could use Gaussian identities to get the marginal distribution
of test outputs $\boldsymbol{g}_{\ast}$.

######  Lemma 9.4 (Marginal Distribution of Test Outputs)

Follow from Section 8.3 (p. 8.3), assume zero-mean prior
$\boldsymbol{\beta}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}_{0})$
and define
$g(\bm{x})=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\beta}+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\epsilon}$
where
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\epsilon\sim\mathcal{N}(0,\sigma^{2})}$.
Given observed training inputs $\bm{X}$, training outputs $\bm{y}$ and test
inputs $\bm{x}_{\ast}$, the marginal distribution of test outputs
$\boldsymbol{g}_{\ast}$ is

$\displaystyle\bm{g}_{\ast}|\bm{X}_{\ast}\bm{X},\bm{y}\sim\mathcal{N}($
$\displaystyle
K(\bm{X}_{\ast},\bm{X})(K(\bm{X},\bm{X})+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\sigma^{2}\bm{I}})^{-1}\bm{y},$
$\displaystyle
K(\bm{X}_{\ast},\bm{x}_{\ast})-K(\bm{X}_{\ast},\bm{X})(K(\bm{X},\bm{X})+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\sigma^{2}\bm{I}})^{-1}K(\bm{X},\bm{X}_{\ast})).$

The marginal distribution of test outputs $\boldsymbol{g}_{\ast}$ agrees with
the form of Equation (9.3) except the output is a vector now (i.e.,
$g_{\ast}\rightarrow\bm{g}_{\ast}$).

####  9.4.3 Further Extension, Generalized Gaussian Process

From the idea of the generalized linear model in Section 6.1 (p. 6.1), we
could set a different noise variance value for each observation. Suppose the
noise covariance matrix now is $\sigma^{2}\boldsymbol{\Lambda}$ where
$\boldsymbol{\Lambda}$ is a diagonal matrix. Then the marginal distribution of
test outputs $\bm{g}_{\ast}$ is

$\displaystyle\bm{g}_{\ast}|\bm{X}_{\ast}\bm{X},\bm{y}\sim\mathcal{N}($
$\displaystyle
K(\bm{X}_{\ast},\bm{X})(K(\bm{X},\bm{X})+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\sigma^{2}\boldsymbol{\Lambda}})^{-1}\bm{y},$
$\displaystyle
K(\bm{X}_{\ast},\bm{x}_{\ast})-K(\bm{X}_{\ast},\bm{X})(K(\bm{X},\bm{X})+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\sigma^{2}\boldsymbol{\Lambda}})^{-1}K(\bm{X},\bm{X}_{\ast})).$

Apparently, the noise covariance should not be dependent on the input matrix
$\bm{X}$. Otherwise, we will not be able to use the kernel trick.

(a) Use half data set

(b) Use all data set

Figure 9.1: Comparison of using half of the data and using all of the data.

(a) Realization 1

(b) Realization 2

(c) Realization 3

(d) Realization 4

Figure 9.2: Random realizations.

(a) Repeat the data set once

(b) Repeat the data set twice

(c) Repeat the data set three times

(d) Repeat the data set four times

Figure 9.3: Repeat data set to have larger data set number.

###### Example 9.4.1

For a data set with input being the house area and output being the rent of
the house. We use cross-validation to select the parameters in the Gaussian
kernel.

In Figure 9.1, the red dots are the training inputs, the blue line is the MAP
estimator via GP regressor and the shaded area are the ones with a confidence
interval of 95%.

In Figure 9.1(b), we use all the data. While in Figure 9.1(a), we only use the
data with house areas smaller than 100 $m^{2}$. We can see the estimator does
not predict well when the house area is larger than 100 $m^{2}$.

In Figure 9.2, we use all the data to predict and we draw four random
realizations from the posterior as shown in orange lines. Thus we can have a
further taste of what the distribution over functions means.

In Figure 9.3, we repeat the data set once, twice, three times, and four times
to see if the number of data sets matter. The yellow line is the predictor by
OLS. We find that the estimate from OLS does not change because it is exactly
the same data set. To see this, suppose the OLS estimator for data without
repeating is given by

$\hat{\boldsymbol{\beta}}_{1}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top}\bm{y}.$

The OLS estimator for data with repeating twice is given by

$\displaystyle\hat{\boldsymbol{\beta}}_{2}$
$\displaystyle=\left(\begin{bmatrix}\bm{X}\\\
\bm{X}\end{bmatrix}^{\top}\begin{bmatrix}\bm{X}\\\
\bm{X}\end{bmatrix}\right)^{-1}\begin{bmatrix}\bm{X}\\\
\bm{X}\end{bmatrix}^{\top}\begin{bmatrix}\bm{y}\\\ \bm{y}\end{bmatrix}$
$\displaystyle=\frac{1}{2}(\bm{X}^{\top}\bm{X})^{-1}\cdot
2\bm{X}^{\top}\bm{y}=\hat{\boldsymbol{\beta}}_{1}.$

which is exactly the same as $\hat{\boldsymbol{\beta}}_{1}$. Similarly, we can
also show that the OLS estimator for data with repeating three times and four
times are the same as well.

However, when we repeat the data set for the GP regressor, the estimate from
the GP regressor has a smaller variance when the data set number increases. In
Section 8.2 (p. 8.2), we mentioned that the Bayesian model considers prior
information on the parameters in the model making it particularly useful to
regularize regression problems where data information is limited. So the data
set number matters in the Bayesian approach where we also gave an example in
Section 8.2 (p. 8.2). This exact example shows the Bayesian thoughts behind
the GP.

###  9.5 Limitations of Gaussian Process from Non Zero-Mean Prior*

When the prior mean $\boldsymbol{\beta}_{0}$ is not zero and the model is
under Gaussian noise
$\boldsymbol{\epsilon}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\bm{I})$, we
transform the predictive prediction into

$\displaystyle g_{\ast}|\bm{x}_{\ast}\bm{X},\bm{y},\sigma^{2}\sim\mathcal{N}($
$\displaystyle\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top}(\bm{K}+\sigma^{2}\bm{I})^{-1}\bm{y}+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{\phi}_{\ast}^{\top}(\frac{1}{\sigma^{2}}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top}\boldsymbol{\Phi}+\bm{I})^{-1}\boldsymbol{\beta}_{0}},$
$\displaystyle\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}_{\ast}-\boldsymbol{\phi}_{\ast}^{\top}\boldsymbol{\Sigma}_{0}\boldsymbol{\Phi}^{\top}(\bm{K}+\sigma^{2}\bm{I})^{-1}\boldsymbol{\Phi}\boldsymbol{\Sigma}_{0}\boldsymbol{\phi}_{\ast}),$

where we use the fact that $(\bm{A}\bm{B})^{-1}=\bm{B}^{-1}\bm{A}^{-1}$.

Following from Section 8.3 (p. 8.3), if we assume a non zero-mean prior
$\boldsymbol{\beta}\sim\mathcal{N}(\boldsymbol{\beta}_{0},\boldsymbol{\Sigma}_{0})$
and define $g(\bm{x})=\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\beta}$.
Then, the mean and covariance for the prior output are:

$\displaystyle\mathrm{E}[g(\bm{x})]$
$\displaystyle=\boldsymbol{\phi}(\bm{x})^{\top}\mathrm{E}[\boldsymbol{\beta}]={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{\phi}(\bm{x})^{\top}\boldsymbol{\beta}_{0}},$
$\displaystyle\mathrm{E}[g(\bm{x})g(\bm{x}^{\prime})]$
$\displaystyle=\boldsymbol{\phi}(\bm{x})^{\top}\mathrm{E}[\boldsymbol{\beta}\boldsymbol{\beta}^{\top}]\boldsymbol{\phi}(\bm{x}^{\prime})=\boldsymbol{\phi}(\bm{x})^{\top}(\boldsymbol{\Sigma}_{0}+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boldsymbol{\beta}_{0}\boldsymbol{\beta}_{0}^{\top}})\boldsymbol{\phi}(\bm{x}^{\prime}),$

Therefore, we could get the joint distribution of the training outputs
$\bm{y}$ and the test outputs $\boldsymbol{g}_{\ast}$

$\left[\begin{matrix}\bm{y}\\\
\boldsymbol{g}_{\ast}\end{matrix}\right]\sim\mathcal{N}\left(\left[\begin{matrix}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi(\bm{x})^{\top}\boldsymbol{\beta}_{0}}\\\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi(\boldsymbol{g}_{\ast})^{\top}\boldsymbol{\beta}_{0}}\end{matrix}\right],\left[\begin{matrix}K(\bm{X},\bm{X})+\sigma^{2}\bm{I},&K(\bm{X},\boldsymbol{X}_{\ast})\\\
K(\boldsymbol{X}_{\ast},\bm{X}),&K(\boldsymbol{X}_{\ast},\boldsymbol{X}_{\ast})\end{matrix}\right]\right).$

By Gaussian identities, the marginal distribution of test outputs is

$\displaystyle\bm{g}_{\ast}|\bm{X}_{\ast}\bm{X},\bm{y}\sim\mathcal{N}($
$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi(\boldsymbol{g}_{\ast})^{\top}\boldsymbol{\beta}_{0}}+K(\bm{X}_{\ast},\bm{X})(K(\bm{X},\bm{X})+\sigma^{2}\bm{I})^{-1}(\bm{y}-{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\phi(\bm{x})^{\top}\boldsymbol{\beta}_{0}}),$
$\displaystyle
K(\bm{X}_{\ast},\bm{x}_{\ast})-K(\bm{X}_{\ast},\bm{X})(K(\bm{X},\bm{X})+\sigma^{2}\bm{I})^{-1}K(\bm{X},\bm{X}_{\ast})).$

where we notice that we cannot find an explicit form in the $z$-dimensional
space for the mean of $g(\bm{x})$ and we cannot apply the kernel trick in this
sense. Thus we lose some flexibility from non zero-mean prior.

##  Chapter 10 Linear Algebra

###  A Dimension of Column Space and Row Space

In this appendix, we prove Lemma 10.1 that the dimension of the column space
of a matrix $\bm{X}\in^{n\times p}$ is equal to the dimension of its row
space, i.e., the row rank and the column rank of a matrix $\bm{X}$ are equal.

######  Lemma 10.1 (Dimension of Column Space and Row Space)

The dimension of the column space of a matrix $\bm{X}\in^{n\times p}$ is equal
to the dimension of its row space, i.e., the row rank and the column rank of a
matrix $\bm{X}$ are equal.

Proof [of Lemma 10.1] We first notice that the null space of $\bm{X}$ is
orthogonal complementary to the row space of $\bm{X}$:
$\mathcal{N}(\bm{X})\bot\mathcal{C}(\bm{X}^{\top})$ (where the row space of
$\bm{X}$ is exactly the column space of $\bm{X}^{\top}$), that is, vectors in
the null space of $\bm{X}$ are orthogonal to vectors in the row space of
$\bm{X}$. To see this, suppose $\bm{X}$ has rows
$\bm{a}_{1}^{\top},\bm{a}_{2}^{\top},\ldots,\bm{a}_{n}^{\top}$ and
$\bm{X}=[\bm{a}_{1}^{\top};\bm{a}_{2}^{\top};\ldots;\bm{a}_{n}^{\top}]$. For
any vector $\boldsymbol{\beta}\in\mathcal{N}(\bm{X})$, we have
$\bm{X}\boldsymbol{\beta}=\boldsymbol{0}$, that is,
$[\bm{a}_{1}^{\top}\boldsymbol{\beta};\bm{a}_{2}^{\top}\boldsymbol{\beta};\ldots;\bm{a}_{n}^{\top}\boldsymbol{\beta}]=\boldsymbol{0}$.
And since the row space of $\bm{X}$ is spanned by
$\bm{a}_{1}^{\top},\bm{a}_{2}^{\top},\ldots,\bm{a}_{n}^{\top}$. Then
$\boldsymbol{\beta}$ is perpendicular to any vectors from
$\mathcal{C}(\bm{X}^{\top})$ which means
$\mathcal{N}(\bm{X})\bot\mathcal{C}(\bm{X}^{\top})$.

Now suppose, the dimension of row space of $\bm{X}$ is $r$. Let
$\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{r}$ be a set of vectors in p and form a
basis for the row space. Then the $r$ vectors
$\bm{X}\bm{r}_{1},\bm{X}\bm{r}_{2},\ldots,\bm{X}\bm{r}_{r}$ are in the column
space of $\bm{X}$, furthermore, they are linearly independent. To see this,
suppose we have a linear combination of the $r$ vectors:
$\beta_{1}\bm{A}\bm{r}_{1}+\beta_{2}\bm{A}\bm{r}_{2}+\ldots+\beta_{r}\bm{A}\bm{r}_{r}=0$,
that is,
$\bm{X}(\beta_{1}\bm{r}_{1}+\beta_{2}\bm{r}_{2}+\ldots+\beta_{r}\bm{r}_{r})=0$
and the vector
$\bm{v}=\beta_{1}\bm{r}_{1}+\beta_{2}\bm{r}_{2}+\ldots+\beta_{r}\bm{r}_{r}$ is
in null space of $\bm{X}$. But since
$\\{\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{r}\\}$ is a basis for the row space
of $\bm{X}$, $\bm{v}$ is thus also in the row space of $\bm{X}$. We have shown
that vectors from null space of $\bm{X}$ is perpendicular to vectors from row
space of $\bm{X}$, thus $\bm{v}^{\top}\bm{v}=0$ and
$\beta_{1}=\beta_{2}=\ldots=\beta_{r}=0$. Then
$\bm{X}\bm{r}_{1},\bm{X}\bm{r}_{2},\ldots,\bm{X}\bm{r}_{r}$ are in the column
space of $\bm{X}$ and they are independent which means the dimension of the
column space of $\bm{X}$ is larger than $r$. This result shows that row rank
of $\bm{X}\leq$ column rank of $\bm{X}$.

If we apply this process again for $\bm{X}^{\top}$. We will have column rank
of $\bm{X}\leq$ row rank of $\bm{X}$. This completes the proof.

Further information can be drawn from this proof is that if
$\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{r}$ is a set of vectors in p that forms
a basis for the row space, then
$\bm{X}\bm{r}_{1},\bm{X}\bm{r}_{2},\ldots,\bm{X}\bm{r}_{r}$ is a basis for the
column space of $\bm{X}$. We formulate this finding into the following lemma.

######  Lemma 10.2 (Column Basis from Row Basis)

For any matrix $\bm{X}\in^{n\times p}$, let
$\\{\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{r}\\}$ be a set of vectors in p which
forms a basis for the row space, then
$\\{\bm{X}\bm{r}_{1},\bm{X}\bm{r}_{2},\ldots,\bm{X}\bm{r}_{r}\\}$ is a basis
for the column space of $\bm{X}$.

###  B Fundamental Theorem of Linear Algebra

Proof [of Theorem 1.4, p. 1.4] Following the proof of Lemma 10.1 in Appendix
A. Let $\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{r}$ be a set of vectors in p that
form a basis for the row space, then
$\bm{X}\bm{r}_{1},\bm{X}\bm{r}_{2},\ldots,\bm{X}\bm{r}_{r}$ is a basis for the
column space of $\bm{X}$. Let $\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{k}\in^{p}$
form a basis for the null space of $\bm{X}$. Following again the proof of
Lemma 10.1, $\mathcal{N}(\bm{X})\bot\mathcal{C}(\bm{X}^{\top})$, thus,
$\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{r}$ are perpendicular to
$\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{k}$. Then,
$\\{\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{r},\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{k}\\}$
is linearly independent in p.

For any vector $\boldsymbol{\beta}\in^{p}$, $\bm{X}\boldsymbol{\beta}$ is in
the column space of $\bm{X}$. Then it can be expressed by a combination of
$\bm{X}\bm{r}_{1},\bm{X}\bm{r}_{2},\ldots,\bm{X}\bm{r}_{r}$:
$\bm{X}\boldsymbol{\beta}=\sum_{i=1}^{r}a_{i}\bm{X}\bm{r}_{i}$ which states
that $\bm{X}(\boldsymbol{\beta}-\sum_{i=1}^{r}a_{i}\bm{r}_{i})=\boldsymbol{0}$
and $\boldsymbol{\beta}-\sum_{i=1}^{r}a_{i}\bm{r}_{i}$ is thus in
$\mathcal{N}(\bm{X})$. Since $\\{\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{k}\\}$
is a basis for the null space of $\bm{X}$,
$\boldsymbol{\beta}-\sum_{i=1}^{r}a_{i}\bm{r}_{i}$ can be expressed by a
combination of $\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{k}$:
$\boldsymbol{\beta}-\sum_{i=1}^{r}a_{i}\bm{r}_{i}=\sum_{j=1}^{k}b_{j}\bm{n}_{j}$,
i.e.,
$\boldsymbol{\beta}=\sum_{i=1}^{r}a_{i}\bm{r}_{i}+\sum_{j=1}^{k}b_{j}\bm{n}_{j}$.
That is, any vector $\boldsymbol{\beta}\in^{p}$ can be expressed by
$\\{\bm{r}_{1},\bm{r}_{2},\ldots,\bm{r}_{r},\bm{n}_{1},\bm{n}_{2},\ldots,\bm{n}_{k}\\}$
and the set forms a basis for p. Thus the dimension sum to $p$: $r+k=p$, i.e.,
$dim(\mathcal{N}(\bm{X}))+dim(\mathcal{C}(\bm{X}^{\top}))=p$. Similarly, we
can prove $dim(\mathcal{N}(\bm{X}^{\top}))+dim(\mathcal{C}(\bm{X}))=n$.

###  C Similar Matrices

######  Definition 10.1 (Similar Matrices)

For any nonsingular matrix $\bm{P}$, the matrices $\bm{A}$ and
$\bm{P}\bm{A}\bm{P}^{-1}$ are called similar matrices.

######  Lemma 10.3 (Eigenvalue and Rank of Similar Matrices)

Any eigenvalue of $\bm{A}$ is also an eigenvalue of $\bm{P}\bm{A}\bm{P}^{-1}$
given any nonsingular matrix $\bm{P}$. The converse is also true that any
eigenvalue of $\bm{P}\bm{A}\bm{P}^{-1}$ is also an eigenvalue of $\bm{A}$.

And also the rank of $\bm{A}$ is equal to the rank of matrix
$\bm{P}\bm{A}\bm{P}^{-1}$ given any nonsingular matrix $\bm{P}$.

Proof [of Lemma 10.3] For any eigenvalue $\lambda$ of $\bm{A}$, we have
$\bm{A}\bm{x}=\lambda\bm{x}$. Then
$\lambda\bm{P}\bm{x}=\bm{P}\bm{A}\bm{P}^{-1}\bm{P}\bm{x}$ such that
$\bm{P}\bm{x}$ is an eigenvector of $\bm{P}\bm{A}\bm{P}^{-1}$ corresponding to
$\lambda$.

Similarly, for any eigenvalue $\lambda$ of $\bm{P}\bm{A}\bm{P}^{-1}$, we have
$\bm{P}\bm{A}\bm{P}^{-1}\bm{x}=\lambda\bm{x}$. Then
$\bm{A}\bm{P}^{-1}\bm{x}=\lambda\bm{P}^{-1}\bm{x}$ such that
$\bm{P}^{-1}\bm{x}$ is an eigenvector of $\bm{A}$ corresponding to $\lambda$.

For the rank of $\bm{P}\bm{A}\bm{P}^{-1}$, we have
$trace(\bm{P}\bm{A}\bm{P}^{-1})=trace(\bm{A}\bm{P}^{-1}\bm{P})=trace(\bm{A})$,
where the first equality comes from the fact that the trace of a product is
invariant under cyclical permutations of the factors:

$\boxed{trace(\bm{A}\bm{B}\bm{C})=trace(\bm{B}\bm{C}\bm{A})=trace(\bm{C}\bm{A}\bm{B})},$

if all $\bm{A}\bm{B}\bm{C}$, $\bm{B}\bm{C}\bm{A}$, and $\bm{C}\bm{A}\bm{B}$
exist.

###  D Cochran’s Theorem

In this appendix, we provide proof for Theorem 4.5 (p. 4.5). To prove the
Cochran’s theorem, we need the following lemma.

######  Lemma 10.4 (Idempotent Decomposition: Rank-Additivity)

For $n\times n$ square matrices $\bm{A}_{1},\bm{A}_{2},\ldots,\bm{A}_{k}$, and
$\bm{A}_{1}+\bm{A}_{2}+\ldots+\bm{A}_{k}=\bm{I}_{n}$, then the following three
conditions are equivalent:

i). $\bm{A}_{i}^{2}=\bm{A}_{i}$, for all $i\in\\{1,2,\ldots,k\\}$, i.e.,
$\bm{A}_{i}$’s are idempotent;

ii). $rank(\bm{A}_{1})+rank(\bm{A}_{2})+\ldots+rank(\bm{A}_{k})=n$;

iii). $\bm{A}_{i}\bm{A}_{j}=\boldsymbol{0}$ for all $i\neq j$ and
$i,j\in\\{1,2,\ldots,k\\}$.

Proof [of Lemma 10.4] From i) to ii), by Lemma 1.19 (p. 1.19), the trace and
rank of any idempotent matrix are the same. Then,

$\sum_{i=1}^{k}rank(\bm{A}_{i})=\sum_{i=1}^{k}trace(\bm{A}_{i})=trace(\bm{I}_{n})=n.$

From ii) to iii), we have the following block Gaussian elimination for a
$(k+1)\times(k+1)$ block matrix (that is, $(k+1)n\times(k+1)n$ matrix) where
$(2,2),(3,3),\ldots,(k+1,k+1)$ blocks are
$\bm{A}_{1},\bm{A}_{2},\ldots,\bm{A}_{k}$ respectively.

$\displaystyle\bm{X}=$ $\displaystyle\begin{bmatrix}\boldsymbol{0}_{n}&&&\\\
&\bm{A}_{1}&&\\\ &&\ddots&\\\ &&&\bm{A}_{k}\\\
\end{bmatrix}\stackrel{{\scriptstyle\bm{E}_{1}\times}}{{\rightarrow}}\begin{bmatrix}\boldsymbol{0}_{n}&\bm{A}_{1}&\ldots&\bm{A}_{k}\\\
&\bm{A}_{1}&\\\ &&\ddots&\\\ &&&\bm{A}_{k}\\\
\end{bmatrix}\stackrel{{\scriptstyle\bm{E}_{2}\times}}{{\rightarrow}}\begin{bmatrix}\bm{I}_{n}&\bm{A}_{1}&\ldots&\bm{A}_{k}\\\
\bm{A}_{1}&\bm{A}_{1}&\\\ \vdots&&\ddots&\\\ \bm{A}_{k}&&&\bm{A}_{k}\\\
\end{bmatrix}\stackrel{{\scriptstyle\bm{E}_{3}\times}}{{\rightarrow}}$
$\displaystyle\begin{bmatrix}\bm{I}_{n}&\bm{A}_{1}&\ldots&\bm{A}_{k}\\\
&\bm{A}_{1}-\bm{A}_{1}^{2}&&-\bm{A}_{1}\bm{A}_{k}\\\ &\vdots&\ddots&\vdots\\\
&-\bm{A}_{k}\bm{A}_{1}&&\bm{A}_{k}-\bm{A}_{k}^{2}\\\
\end{bmatrix}\stackrel{{\scriptstyle\bm{E}_{4}\times}}{{\rightarrow}}\begin{bmatrix}\bm{I}_{n}&&&\\\
&\bm{A}_{1}-\bm{A}_{1}^{2}&&-\bm{A}_{1}\bm{A}_{k}\\\ &\vdots&\ddots&\vdots\\\
&-\bm{A}_{k}\bm{A}_{1}&&\bm{A}_{k}-\bm{A}_{k}^{2}\\\ \end{bmatrix}=\bm{Y},$

where blank entries indicate zeros. And

$\bullet$ $\bm{E}_{1}$ is adding row-2, row-3, $\ldots$, row-(k+1) to row-1;

$\bullet$ $\bm{E}_{2}$ is adding column-2, column-3, $\ldots$, column-(k+1) to
column-1;

$\bullet$ $\bm{E}_{3}$ is subtracting the row-2 by $\bm{A}_{1}\times$(row-1),
subtracting the row-3 by $\bm{A}_{2}\times$(row-1), $\ldots$;

$\bullet$ $\bm{E}_{4}$ is subtracting the column-2 by
$\bm{A}_{1}\times$(column-1), subtracting column-3 by
$\bm{A}_{2}\times$(column-1), $\ldots$.

We notice that elementary operations/transformations will not change the rank
of the matrix. Since $\bm{X}$ is of rank $\sum_{i=1}^{k}r_{k}=n$, and
$\bm{I}_{n}$ in $\bm{Y}$ is of rank-$n$ as well. We must have

$\begin{bmatrix}\bm{A}_{1}-\bm{A}_{1}^{2}&&-\bm{A}_{1}\bm{A}_{k}\\\
\vdots&\ddots&\vdots\\\ -\bm{A}_{k}\bm{A}_{1}&&\bm{A}_{k}-\bm{A}_{k}^{2}\\\
\end{bmatrix}=\boldsymbol{0},$

which implies $\bm{A}_{i}\bm{A}_{j}=\boldsymbol{0}$ for all $i\neq j$ and
$i,j\in\\{1,2,\ldots,k\\}$.

From iii) to i). we have

$\displaystyle\bm{A}_{i}$ $\displaystyle=\bm{A}_{i}\bm{I}_{n}$
$\displaystyle=\bm{A}_{i}(\bm{A}_{1}+\bm{A}_{2}+\ldots+\bm{A}_{k})$
$\displaystyle=\bm{A}_{i}^{2}.$

This completes the proof.

##### Note on nomenclature

We say that the $\bm{A}_{i}$’s are orthogonal if iii) of Lemma 10.4 is
satisfied and the rank is additive if ii) is satisfied.

Now we are ready to prove the Cochran’s Theorem as follows:

Proof [of Theorem 4.5, p. 4.5] From Lemma 10.4, $\bm{A}_{i}$’s are idempotent.

Case 1, If $\bm{y}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\bm{I})$:

By Spectral Theorem 1.16 (p. 1.16) and Lemma 1.17 (p. 1.17, the only possible
eigenvalues of idempotent matrices are 0 and 1), we can rewrite the $q_{i}$ by
$q_{i}=\bm{y}^{\top}\bm{A}_{i}\bm{y}=\bm{y}^{\top}(\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{\top})\bm{y}$,
where $\bm{A}_{i}=\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{\top}$ is the spectral
decomposition of $\bm{A}_{i}$. From the fact that rotations on the normal
distribution do not affect the distribution111Rotations on the Gaussian
distribution do not affect the distribution. That is for any orthogonal matrix
$\bm{Q}$ with $\bm{Q}\bm{Q}^{\top}=\bm{Q}^{\top}\bm{Q}=\bm{I}$, if
$\bm{X}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\bm{I})$, then
$\bm{Q}\bm{X}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\bm{I})$., we can define

$\boldsymbol{\eta}=\bm{Q}^{\top}\bm{y}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\bm{I}).$

Thus,

$q_{i}=\boldsymbol{\eta}^{\top}\boldsymbol{\Lambda}\boldsymbol{\eta}\sim\sigma^{2}\chi^{2}_{\mathrm{rank(\bm{A}_{i})}}\sim\sigma^{2}\chi^{2}_{(r_{i})},$

Case 2, If $\bm{y}\sim\mathcal{N}(\boldsymbol{\mu},\sigma^{2}\bm{I})$, and
$\boldsymbol{\mu}^{\top}\bm{A}_{i}\boldsymbol{\mu}=0$:

Let
$p_{i}=(\bm{y}-\boldsymbol{\mu})^{\top}\bm{A}_{i}(\bm{y}-\boldsymbol{\mu})$,
similarly, we can rewrite the $p_{i}$ by
$p_{i}=(\bm{y}-\boldsymbol{\mu})^{\top}\bm{A}_{i}(\bm{y}-\boldsymbol{\mu})=(\bm{y}-\boldsymbol{\mu})^{\top}(\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{\top})(\bm{y}-\boldsymbol{\mu})$,
where $\bm{A}_{i}=\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{\top}$ is the spectral
decomposition of $\bm{A}_{i}$. From the fact that rotations on the normal
distribution do not effect the distribution, we can define

$\boldsymbol{\eta}=\bm{Q}^{\top}(\bm{y}-\boldsymbol{\mu})\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}\bm{I}).$

Thus,

$p_{i}=\boldsymbol{\eta}^{\top}\boldsymbol{\Lambda}\boldsymbol{\eta}\sim\sigma^{2}\chi^{2}_{\mathrm{rank(\bm{A}_{i})}}\sim\sigma^{2}\chi^{2}_{(r_{i})}.$

Then, we decompose the $p_{i}$ into

$\displaystyle p_{i}$
$\displaystyle=\bm{y}^{\top}\bm{A}_{i}\bm{y}-2\bm{y}^{\top}\bm{A}_{i}\boldsymbol{\mu}+\boldsymbol{\mu}^{\top}\bm{A}_{i}\boldsymbol{\mu}$
$\displaystyle=q_{i}-2\bm{y}^{\top}(\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{\top})\boldsymbol{\mu}+\boldsymbol{\mu}^{\top}(\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{\top})\boldsymbol{\mu}.$

Since we assume
$\boldsymbol{\mu}^{\top}\bm{A}_{i}\boldsymbol{\mu}=\boldsymbol{\mu}^{\top}(\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{\top})\boldsymbol{\mu}=0$,
and $\boldsymbol{\Lambda}$ contains only 1 and 0 on the diagonal. We have
$\boldsymbol{\Lambda}=\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{\top}$. That
is

$\boldsymbol{\mu}^{\top}(\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{\top})\boldsymbol{\mu}=\boldsymbol{\mu}^{\top}(\bm{Q}\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{\top}\bm{Q}^{\top})\boldsymbol{\mu}=||\boldsymbol{\Lambda}^{\top}\bm{Q}^{\top}\boldsymbol{\mu}||^{2}=0,$

which implies
$2\bm{y}^{\top}\bm{A}_{i}\boldsymbol{\mu}=2\bm{y}^{\top}(\bm{Q}\boldsymbol{\Lambda}\bm{Q}^{\top})\boldsymbol{\mu}=0$.
Thus $q_{i}=q_{i}\sim\sigma^{2}\chi^{2}_{(r_{i})}$.

Case 3: From Lemma 10.4, $\bm{A}_{i}\bm{A}_{j}=\boldsymbol{0}$ if $i\neq j$.
Let $\bm{A}_{i}=\bm{Q}_{i}\boldsymbol{\Lambda}_{i}\bm{Q}_{i}^{\top}$,
$\bm{A}_{j}=\bm{Q}_{j}\boldsymbol{\Lambda}_{j}\bm{Q}_{j}^{\top}$ be the
spectral decomposition of $\bm{A}_{i}$ and $\bm{A}_{j}$. Then, we have

$\displaystyle\bm{A}_{i}\bm{A}_{j}$
$\displaystyle=\bm{Q}_{i}\boldsymbol{\Lambda}_{i}\bm{Q}_{i}^{\top}\bm{Q}_{j}\boldsymbol{\Lambda}_{j}\bm{Q}_{j}^{\top}=\boldsymbol{0}$
(10.1) $\displaystyle\bm{Q}_{i}^{\top}\bm{A}_{i}\bm{A}_{j}\bm{Q}_{j}$
$\displaystyle=\boldsymbol{\Lambda}_{i}\bm{Q}_{i}^{\top}\bm{Q}_{j}\boldsymbol{\Lambda}_{j}=\boldsymbol{0}.$

Write out $q_{i}$ and $q_{j}$:

$\displaystyle q_{i}$
$\displaystyle=\bm{y}^{\top}\bm{Q}_{i}\boldsymbol{\Lambda}_{i}\bm{Q}_{i}^{\top}\bm{y}=\bm{y}^{\top}\bm{Q}_{i}\boldsymbol{\Lambda}_{i}\boldsymbol{\Lambda}_{i}^{\top}\bm{Q}_{i}^{\top}\bm{y}$
$\displaystyle q_{j}$
$\displaystyle=\bm{y}^{\top}\bm{Q}_{j}\boldsymbol{\Lambda}_{j}\bm{Q}_{j}^{\top}\bm{y}=\bm{y}^{\top}\bm{Q}_{j}\boldsymbol{\Lambda}_{j}\boldsymbol{\Lambda}_{j}^{\top}\bm{Q}_{j}^{\top}\bm{y}.$

Let $\bm{a}_{i}=\boldsymbol{\Lambda}_{i}^{\top}\bm{Q}_{i}^{\top}\bm{y}$ and
$\bm{a}_{j}=\boldsymbol{\Lambda}_{j}^{\top}\bm{Q}_{j}^{\top}\bm{y}$, we have

$\displaystyle\mathrm{Cov}[\bm{a}_{i},\bm{a}_{j}]$
$\displaystyle=\boldsymbol{\Lambda}_{i}^{\top}\bm{Q}_{i}^{\top}\mathrm{Cov}[\bm{y},\bm{y}]\bm{Q}_{j}\boldsymbol{\Lambda}_{j}=\sigma^{2}\boldsymbol{\Lambda}_{i}^{\top}\bm{Q}_{i}^{\top}\bm{Q}_{j}\boldsymbol{\Lambda}_{j}=\boldsymbol{0},$

where the last equality is from Equation (10.1). This implies
$\mathrm{Cov}[q_{i},q_{j}]=0$ since $q_{i}=\bm{a}_{i}^{\top}\bm{a}_{i}$ and
$q_{j}=\bm{a}_{j}^{\top}\bm{a}_{j}$. And we complete the proof.
Another proof is provided in Gut (2009b). But the author did not provide an
inductive case for $k>2$. Interesting readers can refer to it.

##  Chapter 11 Matrix Decomposition

###  A CR Decomposition

The CR decomposition is proposed in Strang (2021); Strang and Moler (2022). We
firstly give the result and we will discuss the existence and the origin of
this decomposition in the following sections.

######  Theorem 11.1 (CR Decomposition)

Any rank-$r$ matrix $\bm{A}\in^{n\times p}$ can be factored as

$\underset{n\times p}{\bm{A}}=\underset{n\times
r}{\bm{C}}\,\,\,\,\,\,\,\,\underset{r\times p}{\bm{R}}$

where $\bm{C}$ is the first $r$ linearly independent columns of $\bm{A}$, and
$\bm{R}$ is an $r\times p$ matrix to reconstruct the columns of $\bm{A}$ from
columns of $\bm{C}$. In particular, $\bm{R}$ is the row reduced echelon form
(RREF) of $\bm{A}$ without the zero rows.

The storage for the decomposition is then reduced or potentially increased
from $np$ to $r(n+p)$.

####  A.1 Existence of the CR Decomposition

Since matrix $\bm{A}$ is of rank $r$, there are some $r$ linearly independent
columns in $\bm{A}$. We then choose linearly independent columns from $\bm{A}$
and put them into $\bm{C}$:

Find $r$ linearly Independent Columns From $\bm{A}$ 1\. If column 1 of
$\bm{A}$ is not zero, put it into the column of $\bm{C}$; 2\. If column 2 of
$\bm{A}$ is not a multiple of column 1, put it into the column of $\bm{C}$;
3\. If column 3 of $\bm{A}$ is not a combination of columns 1 and 2, put it
into the column of $\bm{C}$; 4\. Continue this process until we find $r$
linearly independent columns (or all the linearly independent columns if we do
not know the rank $r$ beforehand).

When we have the $r$ linearly independent columns from $\bm{A}$, we can prove
the existence of CR decomposition by the column space view of matrix
multiplication.

##### Column space view of matrix multiplication

A multiplication of two matrices $\bm{D}\in^{n\times k},\bm{E}\in^{k\times p}$
is
$\bm{A}=\bm{D}\bm{E}=\bm{D}[\bm{e}_{1},\bm{e}_{2},\ldots,\bm{e}_{p}]=[\bm{D}\bm{e}_{1},\bm{D}\bm{e}_{2},\ldots,\bm{D}\bm{e}_{p}]$,
i.e., each column of $\bm{A}$ is a combination of columns from $\bm{D}$.

Proof [of Theorem 11.1] As the rank of matrix $\bm{A}$ is $r$ and $\bm{C}$
contains $r$ linearly independent columns from $\bm{A}$, the column space of
$\bm{C}$ is equivalent to the column space of $\bm{A}$. If we take any other
column $\bm{a}_{i}$ of $\bm{A}$, $\bm{a}_{i}$ can be represented as a linear
combination of the columns of $\bm{C}$, i.e., there exists a vector
$\bm{r}_{i}$ such that $\bm{a}_{i}=\bm{C}\bm{r}_{i}$, $\forall
i\in\\{1,2,\ldots,p\\}$. Put these $\bm{r}_{i}$’s into the columns of matrix
$\bm{R}$, we obtain

$\bm{A}=[\bm{a}_{1},\bm{a}_{2},\ldots,\bm{a}_{p}]=[\bm{C}\bm{r}_{1},\bm{C}\bm{r}_{2},\ldots,\bm{C}\bm{r}_{p}]=\bm{C}\bm{R},$

from which the result follows.

####  A.2 Reduced Row Echelon Form (RREF)

We write out the Gaussian elimination for a $4\times 4$ square matrix, where
$\boxtimes$ represents a value that is not necessarily zero, and boldface
indicates the value has just been changed:

Gaussian Elimination for a Square Matrix
$\mathop{{}\begin{bmatrix}\boxtimes&\boxtimes&\boxtimes&\boxtimes\\\
\boxtimes&\boxtimes&\boxtimes&\boxtimes\\\
\boxtimes&\boxtimes&\boxtimes&\boxtimes\\\
\boxtimes&\boxtimes&\boxtimes&\boxtimes\end{bmatrix}}_{\textstyle\mathstrut\bm{A}}\stackrel{{\scriptstyle\bm{E}_{1}}}{{\longrightarrow}}\mathop{{}\begin{bmatrix}\boxtimes&\boxtimes&\boxtimes&\boxtimes\\\
\bm{0}&\bm{0}&\bm{\boxtimes}&\bm{\boxtimes}\\\
\bm{0}&\bm{\boxtimes}&\bm{\boxtimes}&\bm{\boxtimes}\\\
\bm{0}&\bm{\boxtimes}&\bm{\boxtimes}&\bm{\boxtimes}\end{bmatrix}}_{\textstyle\mathstrut\bm{E}_{1}\bm{A}}\stackrel{{\scriptstyle\bm{P}_{1}}}{{\longrightarrow}}\mathop{{}\begin{bmatrix}\boxtimes&\boxtimes&\boxtimes&\boxtimes\\\
\bm{0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{\boxtimes}}&\bm{\boxtimes}&\bm{\boxtimes}\\\
\bm{0}&\bm{0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{\boxtimes}}&\bm{\boxtimes}\\\
0&\boxtimes&\boxtimes&\boxtimes\end{bmatrix}}_{\textstyle\mathstrut\bm{P}_{1}\bm{E}_{1}\bm{A}}\stackrel{{\scriptstyle\bm{E}_{2}}}{{\longrightarrow}}\mathop{{}\begin{bmatrix}\boxtimes&\boxtimes&\boxtimes&\boxtimes\\\
0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boxtimes}&\boxtimes&\boxtimes\\\
0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\boxtimes}&\boxtimes\\\
0&\bm{0}&\bm{0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{\boxtimes}}\end{bmatrix}}_{\textstyle\mathstrut\bm{E}_{2}\bm{P}_{1}\bm{E}_{1}\bm{A}}.$

Furthermore, the Gaussian elimination can also be applied to a rectangular
matrix, we give an example for a $4\times 5$ matrix as follows:

Gaussian Elimination for a Rectangular Matrix
$\mathop{{}\begin{bmatrix}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2}&\boxtimes&10&9&\boxtimes\\\
\boxtimes&\boxtimes&\boxtimes&\boxtimes&\boxtimes\\\
\boxtimes&\boxtimes&\boxtimes&\boxtimes&\boxtimes\\\
\boxtimes&\boxtimes&\boxtimes&\boxtimes&\boxtimes\\\
\end{bmatrix}}_{\textstyle\mathstrut\bm{A}}\stackrel{{\scriptstyle\bm{E}_{1}}}{{\longrightarrow}}\mathop{{}\begin{bmatrix}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2}&\boxtimes&10&9&\boxtimes\\\
\bm{0}&\bm{0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{5}}&\bm{6}&\bm{\boxtimes}\\\
\bm{0}&\bm{0}&\bm{2}&\bm{\boxtimes}&\bm{\boxtimes}\\\
\bm{0}&\bm{0}&\bm{\boxtimes}&\bm{\boxtimes}&\bm{\boxtimes}\\\
\end{bmatrix}}_{\textstyle\mathstrut\bm{E}_{1}\bm{A}}\stackrel{{\scriptstyle\bm{E}_{2}}}{{\longrightarrow}}\mathop{{}\begin{bmatrix}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2}&\boxtimes&10&9&\boxtimes\\\
0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}5}&6&\boxtimes\\\
0&0&\bm{0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{3}}&\bm{\boxtimes}\\\
0&0&\bm{0}&\bm{0}&\bm{0}\\\
\end{bmatrix}}_{\textstyle\mathstrut\bm{E}_{2}\bm{E}_{1}\bm{A}},$

where the blue-colored numbers are pivots and we call the last matrix above
row echelon form. Note that we get the 4-th row as a zero row in this specific
example. Going further, we subtract each row by a multiple of the next row to
make the entries above the pivots to be zero:

Reduced Row Echelon Form: Get Zero Above Pivots
$\mathop{{}\begin{bmatrix}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2}&\boxtimes&10&9&\boxtimes\\\
0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}5}&6&\boxtimes\\\
0&0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}3}&\boxtimes\\\
0&0&0&0&0\\\
\end{bmatrix}}_{\textstyle\mathstrut\bm{E}_{2}\bm{E}_{1}\bm{A}}\stackrel{{\scriptstyle\bm{E}_{3}}}{{\longrightarrow}}\mathop{{}\begin{bmatrix}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2}&\boxtimes&\bm{0}&\bm{-3}&\bm{\boxtimes}\\\
0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}5}&6&\boxtimes\\\
0&0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}3}&\boxtimes\\\
0&0&0&0&0\\\
\end{bmatrix}}_{\textstyle\mathstrut\bm{E}_{3}\bm{E}_{2}\bm{E}_{1}\bm{A}}\stackrel{{\scriptstyle\bm{E}_{4}}}{{\longrightarrow}}\mathop{{}\begin{bmatrix}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2}&\boxtimes&0&\bm{0}&\bm{\boxtimes}\\\
0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}5}&\bm{0}&\bm{\boxtimes}\\\
0&0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}3}&\boxtimes\\\
0&0&0&0&0\\\
\end{bmatrix}}_{\textstyle\mathstrut\bm{E}_{4}\bm{E}_{3}\bm{E}_{2}\bm{E}_{1}\bm{A}},$

where $\bm{E}_{3}$ subtracts 2 times the $2$-nd row from the $1$-st row, and
$\bm{E}_{4}$ adds the $3$-rd row to the $1$-st row and subtracts 2 times the
$3$-rd row from the $2$-nd row. Finally, we get the full row reduced echelon
form by making the pivots to be 1:

Reduced Row Echelon Form: Make The Pivots To Be 1
$\mathop{{}\begin{bmatrix}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}2}&\boxtimes&0&0&\boxtimes\\\
0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}5}&0&\boxtimes\\\
0&0&0&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}3}&\boxtimes\\\
0&0&0&0&0\\\
\end{bmatrix}}_{\textstyle\mathstrut\bm{E}_{4}\bm{E}_{3}\bm{E}_{2}\bm{E}_{1}\bm{A}}\stackrel{{\scriptstyle\bm{E}_{5}}}{{\longrightarrow}}\mathop{{}\begin{bmatrix}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{1}}&\bm{\boxtimes}&\bm{0}&\bm{0}&\bm{\boxtimes}\\\
\bm{0}&\bm{0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{1}}&\bm{0}&\bm{\boxtimes}\\\
\bm{0}&\bm{0}&\bm{0}&{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\bm{1}}&\bm{\boxtimes}\\\
0&0&0&0&0\\\
\end{bmatrix}}_{\textstyle\mathstrut\bm{E}_{5}\bm{E}_{4}\bm{E}_{3}\bm{E}_{2}\bm{E}_{1}\bm{A}},$

where $\bm{E}_{5}$ makes the pivots to be 1. We call this final matrix the
reduced row echelon form of $\bm{A}$ where it has 1’s as pivots and zeros
above the pivots.

######  Lemma 11.2 (Rank and Pivots)

The rank of $\bm{A}$ is equal to the number of pivots.

######  Lemma 11.3 (RREF in CR)

The reduced row echelon form of the matrix $\bm{A}$ without zero rows is the
matrix $\bm{R}$ in the CR decomposition.

Proof [Informal Proof of Lemma 11.3 and Lemma 11.2] Following the steps in the
Gaussian elimination, the number of pivot elements indicates the number of
linearly independent columns in matrix $\bm{A}$, which is on the other hand,
exactly the rank of the matrix.

Following the example above, we have

$\bm{E}_{5}\bm{E}_{4}\bm{E}_{3}\bm{E}_{2}\bm{E}_{1}\bm{A}=\bm{R}\quad\longrightarrow\quad\bm{A}=(\bm{E}_{5}\bm{E}_{4}\bm{E}_{3}\bm{E}_{2}\bm{E}_{1})^{-1}\bm{R}.$

We notice that the columns 1, 3, 4 of $\bm{R}$ have only one element 1 which
means we could construct a matrix $\bm{C}$ (exactly the same column matrix in
CR decomposition) whose first 3 columns are the same as the 1-st, 3-rd, 4-th
columns of matrix $\bm{A}$ respectively, i.e.,
$\bm{C}=[\bm{a}_{1},\bm{a}_{3},\bm{a}_{4}]$. Furthermore, since the last row
of $\bm{R}$ is a zero vector, the 4-th column of $\bm{C}$ does not account for
any computation so that we can just ignore the 4-th column of $\bm{C}$ and the
4-th row of $\bm{R}$. And this $\bm{C}$ is the only matrix that can
reconstruct the 1-st, 3-rd, 4-th columns of $\bm{A}$ as the pivots of $\bm{R}$
are all 1. We obtain

$\bm{A}=\bm{C}\bm{R}.$

This completes the proof.
In short, we first compute the reduced row echelon form of the matrix $\bm{A}$
by $rref(\bm{A})$, Then $\bm{C}$ is obtained by removing from $\bm{A}$ all the
non-pivot columns (this can be determined by looking for columns in
$rref(\bm{A})$ which do not contain a pivot). And $\bm{R}$ is obtained by
eliminating zero rows of $rref(\bm{A})$. And this is actually a special case
of rank decomposition of matrix $\bm{A}$. However, CR decomposition is so
special that it involves the reduced row echelon form so we introduce it here
particularly.

$\bm{R}$ has a remarkable form whose $r$ columns containing the pivots form an
$r\times r$ identity matrix. Note again that we can just remove the zero rows
from the row reduced echelon form to obtain this matrix $\bm{R}$. In Strang
(2021), the authors give a specific notation for the row reduced echelon form
without removing the zero rows as $\bm{R}_{0}$:

$\bm{R}_{0}=rref(\bm{A})=\begin{bmatrix}\bm{R}\\\
\boldsymbol{0}\end{bmatrix}=\begin{bmatrix}\bm{I}_{r}&\bm{F}\\\
\boldsymbol{0}&\boldsymbol{0}\end{bmatrix}\bm{P},$

where the $p\times p$ permutation matrix $\bm{P}$ puts the columns of $r\times
r$ identity matrix $\bm{I}_{r}$ into the correct positions, matching the first
$r$ linearly independent columns of the original matrix $\bm{A}$.

The important theorem in linear algebra that the row rank equals the column
rank of any matrix can be proved by the UTV framework (Theorem 10.1, p. 10.1,
see Lu (2021a) for using UTV to prove the theorem). The CR decomposition also
reveals the theorem.

Proof [of Theorem 10.1, p. 10.1, A Second Way] For CR decomposition of matrix
$\bm{A}=\bm{C}\bm{R}$, we have $\bm{R}=[\bm{I}_{r},\bm{F}]\bm{P}$, where
$\bm{P}$ is a $p\times p$ permutation to put the columns of the $r\times r$
identity matrix $\bm{I}_{r}$ into the correct positions as shown above. It can
be easily verified that the $r$ rows of $\bm{R}$ are linearly independent
because of the submatrix $\bm{I}_{r}$ (since $\bm{I}_{r}$ is nonsingular) such
that the row rank of $\bm{R}$ is $r$.

Firstly, from the definition of the CR decomposition, the $r$ columns of
$\bm{C}$ are from $r$ linearly independent columns of $\bm{A}$, and the column
rank of $\bm{A}$ is $r$. Further,

$\bullet$ Since $\bm{A}=\bm{C}\bm{R}$, all rows of $\bm{A}$ are combinations
of the rows of $\bm{R}$. That is, the row rank of $\bm{A}$ is no larger than
the row rank of $\bm{R}$;

$\bullet$ From $\bm{A}=\bm{C}\bm{R}$, we also have
$(\bm{C}^{\top}\bm{C})^{-1}\bm{C}^{\top}\bm{C}\bm{R}=(\bm{C}^{\top}\bm{C})^{-1}\bm{C}^{\top}\bm{A}$,
that is $\bm{R}=(\bm{C}^{\top}\bm{C})^{-1}\bm{C}^{\top}\bm{A}$.
$\bm{C}^{\top}\bm{C}$ is nonsingular since it has full column rank $r$. Then
all rows of $\bm{R}$ are also combinations of the rows of $\bm{A}$. That is,
the row rank of $\bm{R}$ is no larger than the row rank of $\bm{A}$;

$\bullet$ By “sandwiching”, the row rank of $\bm{A}$ is equal to the row rank
of $\bm{R}$ which is $r$.

Therefore, both the row rank and column rank of $\bm{A}$ are equal to $r$ from
which the result follows.

In the above proof, we apply CR decomposition to show that the row rank of a
matrix is equal to its column rank. Moreover, we will also discuss the special
form of pseudo-inverse from CR decomposition in Appendix A (p. A).

####  A.3 Rank Decomposition

We previously mentioned that the CR decomposition is a special case of rank
decomposition. Formally, we prove the existence of the rank decomposition
rigorously in the following theorem.

######  Theorem 11.4 (Rank Decomposition)

Any rank-$r$ matrix $\bm{A}\in^{n\times p}$ can be factored as

$\underset{n\times p}{\bm{A}}=\underset{n\times
r}{\bm{D}}\,\,\,\,\,\,\,\,\underset{r\times p}{\bm{F}},$

where $\bm{D}\in^{n\times r}$ has rank $r$, and $\bm{F}\in^{r\times p}$ also
has rank $r$, i.e., $\bm{D},\bm{F}$ have full rank $r$.

The storage for the decomposition is then reduced or potentially increased
from $np$ to $r(n+p)$.

Proof [of Theorem 11.4] By ULV decomposition in Theorem 1.8 (p. 1.8), we can
decompose $\bm{A}$ by

$\bm{A}=\bm{U}\begin{bmatrix}\bm{L}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}\end{bmatrix}\bm{V}.$

Let $\bm{U}_{0}=\bm{U}_{:,1:r}$ and $\bm{V}_{0}=\bm{V}_{1:r,:}$, i.e.,
$\bm{U}_{0}$ contains only the first $r$ columns of $\bm{U}$, and $\bm{V}_{0}$
contains only the first $r$ rows of $\bm{V}$. Then, we still have
$\bm{A}=\bm{U}_{0}\bm{L}\bm{V}_{0}$ where $\bm{U}_{0}\in^{n\times r}$ and
$\bm{V}_{0}\in^{r\times p}$. This is also known as the reduced ULV
decomposition as shown in Figure 1.6 (p. 1.6). Let {$\bm{D}=\bm{U}_{0}\bm{L}$
and $\bm{F}=\bm{V}_{0}$}, or {$\bm{D}=\bm{U}_{0}$ and
$\bm{F}=\bm{L}\bm{V}_{0}$}, we find such rank decompositions.
The rank decomposition is not unique. Even by elementary transformations, we
have

$\bm{A}=\bm{E}_{1}\begin{bmatrix}\bm{Z}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}\end{bmatrix}\bm{E}_{2},$

where $\bm{E}_{1}\in^{n\times n},\bm{E}_{2}\in^{p\times p}$ represent
elementary row and column operations, $\bm{Z}\in^{r\times r}$. The
transformation is rather general, and there are dozens of these
$\bm{E}_{1},\bm{E}_{2},\bm{Z}$. By similar construction on this decomposition
as shown in the above proof, we can recover another rank decomposition.

Analogously, we can find such $\bm{D},\bm{F}$ by SVD, URV, CR, CUR, and many
other decompositional algorithms. However, we may connect the different rank
decompositions by the following lemma.

######  Lemma 11.5 (Connection Between Rank Decompositions)

For any two rank decompositions of
$\bm{A}=\bm{D}_{1}\bm{F}_{1}=\bm{D}_{2}\bm{F}_{2}$, there exists a nonsingular
matrix $\bm{P}$ such that

$\bm{D}_{1}=\bm{D}_{2}\bm{P}\qquad\text{and}\qquad\bm{F}_{1}=\bm{P}^{-1}\bm{F}_{2}.$

Proof [of Lemma 11.5] Since $\bm{D}_{1}\bm{F}_{1}=\bm{D}_{2}\bm{F}_{2}$, we
have
$\bm{D}_{1}\bm{F}_{1}\bm{F}_{1}^{\top}=\bm{D}_{2}\bm{F}_{2}\bm{F}_{1}^{\top}$.
It is trivial that $rank(\bm{F}_{1}\bm{F}_{1}^{\top})=rank(\bm{F}_{1})=r$ such
that $\bm{F}_{1}\bm{F}_{1}^{\top}$ is a square matrix with full rank and thus
is nonsingular. This implies
$\bm{D}_{1}=\bm{D}_{2}\bm{F}_{2}\bm{F}_{1}^{\top}(\bm{F}_{1}\bm{F}_{1}^{\top})^{-1}$.
Let $\bm{P}=\bm{F}_{2}\bm{F}_{1}^{\top}(\bm{F}_{1}\bm{F}_{1}^{\top})^{-1}$, we
have $\bm{D}_{1}=\bm{D}_{2}\bm{P}$ and $\bm{F}_{1}=\bm{P}^{-1}\bm{F}_{2}$.

###  B QR Decomposition

In this appendix, we provide the proof for Theorem 1.5 (p. 1.5), the existence
of QR decomposition. In many applications, we are interested in the column
space of a matrix
$\bm{X}=[\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{p}]\in^{n\times p}$. The
successive spaces spanned by the columns $\bm{x}_{1},\bm{x}_{2},\ldots$ of
$\bm{X}$ are

$\mathcal{C}([\bm{x}_{1}])\,\,\,\,\subseteq\,\,\,\,\mathcal{C}([\bm{x}_{1},\bm{x}_{2}])\,\,\,\,\subseteq\,\,\,\,\mathcal{C}([\bm{x}_{1},\bm{x}_{2},\bm{x}_{3}])\,\,\,\,\subseteq\,\,\,\,\ldots,$

where $\mathcal{C}([\ldots])$ is the subspace spanned by the vectors included
in the brackets. The idea of QR decomposition is the construction of a
sequence of orthonormal vectors $\bm{q}_{1},\bm{q}_{2},\ldots$ that span the
same successive subspaces. That is

$\mathcal{C}([\bm{q}_{1}])=\mathcal{C}([\bm{x}_{1}]),\qquad\mathcal{C}([\bm{q}_{1},\bm{q}_{2}])=\mathcal{C}([\bm{x}_{1},\bm{x}_{2}]),\qquad\mathcal{C}([\bm{q}_{1},\bm{q}_{2},\bm{q}_{3}])=\mathcal{C}([\bm{x}_{1},\bm{x}_{2},\bm{x}_{3}]),\ldots.$

To achieve this, we first discuss how to project a vector onto another vector,
based on which the Gram-Schmidt process is employed iteratively.

####  B.1 Project a Vector Onto Another Vector

Project a vector $\bm{a}$ to a vector $\bm{b}$ is to find the vector closest
to $\bm{a}$ on the line of $\bm{b}$. The projection vector $\hat{\bm{a}}$ is
some multiple of $\bm{b}$. Let $\hat{\bm{a}}=\hat{x}\bm{b}$ and
$\bm{a}-\hat{\bm{a}}$ is perpendicular to $\bm{b}$ as shown in Figure 11.1(a).
We then get the following result:

Project Vector $\bm{a}$ Onto Vector $\bm{b}$ $\bm{a}-\hat{\bm{a}}$ is
perpendicular to $\bm{b}$, so $(\bm{a}-\hat{x}\bm{b})^{\top}\bm{b}=0$:
$\hat{x}$ = $\frac{\bm{a}^{\top}\bm{b}}{\bm{b}^{\top}\bm{b}}$ and
$\hat{\bm{a}}=\frac{\bm{a}^{\top}\bm{b}}{\bm{b}^{\top}\bm{b}}\bm{b}=\frac{\bm{b}\bm{b}^{\top}}{\bm{b}^{\top}\bm{b}}\bm{a}$.

(a) Project onto a line

(b) Project onto a space

Figure 11.1: Project a vector onto a line and onto a space.

####  B.2 Project a Vector Onto a Plane

Project a vector $\bm{a}$ to a space spanned by
$\bm{b}_{1},\bm{b}_{2},\ldots,\bm{b}_{n}$ is to find the vector closest to
$\bm{a}$ on the column space of $[\bm{b}_{1},\bm{b}_{2},\ldots,\bm{b}_{n}]$.
The projection vector $\widehat{\bm{a}}$ is a combination of
$\bm{b}_{1},\bm{b}_{2},\ldots,\bm{b}_{n}$:
$\widehat{\bm{a}}=\widehat{x}_{1}\bm{b}_{1}+\widehat{x}_{2}\bm{b}_{2}+\ldots+\widehat{x}_{n}\bm{b}_{n}$.
This is actually a least squares problem. To find the projection, we just
solve the normal equation
$\bm{B}^{\top}\bm{B}\widehat{\bm{x}}=\bm{B}^{\top}\bm{a}$ where
$\bm{B}=[\bm{b}_{1},\bm{b}_{2},\ldots,\bm{b}_{n}]$ and
$\widehat{\bm{x}}=[\widehat{x}_{1},\widehat{x}_{2},\ldots,\widehat{x}_{n}]$.
For each vector $\bm{b}_{i}$, the projection of $\bm{a}$ in the direction of
$\bm{b}_{i}$ can be analogously obtained by

$\widehat{\bm{a}}_{i}=\frac{\bm{b}_{i}\bm{b}_{i}^{\top}}{\bm{b}_{i}^{\top}\bm{b}_{i}}\bm{a},\,\,\,\,\,\,\,\,\forall
i\in\\{1,2,\ldots,n\\}.$

Let $\widehat{\bm{a}}=\sum_{i=1}^{n}\widehat{\bm{a}_{i}}$, this results in

$\bm{a}^{\perp}=(\bm{a}-\widehat{\bm{a}})\perp\mathcal{C}(\bm{B}),$

i.e., $(\bm{a}-\widehat{\bm{a}})$ is perpendicular to the column space of
$\bm{B}=[\bm{b}_{1},\bm{b}_{2},\ldots,\bm{b}_{n}]$ as shown in Figure 11.1(b).

####  B.3 Existence of the QR Decomposition via the Gram-Schmidt Process

We then develop the Gram-Schmidt process based on the projection of a vector.
Given three independent vectors $\bm{a}_{1},\bm{a}_{2},\bm{a}_{3}$ and the
space spanned by the three vectors
$\mathcal{C}{([\bm{a}_{1},\bm{a}_{2},\bm{a}_{3}])}$, i.e., the column space of
matrix $[\bm{a}_{1},\bm{a}_{2},\bm{a}_{3}]$. We intend to construct three
orthogonal vectors $\bm{b}_{1},\bm{b}_{2},\bm{b}_{3}$ in which case
$\mathcal{C}{([\bm{b}_{1},\bm{b}_{2},\bm{b}_{3}])}$ =
$\mathcal{C}{([\bm{a}_{1},\bm{a}_{2},\bm{a}_{3}])}$. Then we divide the
orthogonal vectors by their length to normalize. This process produces three
orthonormal vectors $\bm{q}_{1}=\frac{\bm{b}_{1}}{||\bm{b}_{1}||}$,
$\bm{q}_{2}=\frac{\bm{b}_{2}}{||\bm{b}_{2}||}$,
$\bm{q}_{2}=\frac{\bm{b}_{2}}{||\bm{b}_{2}||}$.

For the first vector, we choose $\bm{b}_{1}=\bm{a}_{1}$ directly. The second
vector $\bm{b}_{1}$ must be perpendicular to the first one. This is actually
the vector $\bm{a}_{2}$ subtracting its projection along $\bm{b}_{1}$:

$\displaystyle\bm{b}_{2}$
$\displaystyle=\bm{a}_{2}-\frac{\bm{b}_{1}\bm{b}_{1}^{\top}}{\bm{b}_{1}^{\top}\bm{b}_{1}}\bm{a}_{2}=(\bm{I}-\frac{\bm{b}_{1}\bm{b}_{1}^{\top}}{\bm{b}_{1}^{\top}\bm{b}_{1}})\bm{a}_{2}\qquad$
$\displaystyle(\text{Projection view})$
$\displaystyle=\bm{a}_{2}-\underbrace{\frac{\bm{b}_{1}^{\top}\bm{a}_{2}}{\bm{b}_{1}^{\top}\bm{b}_{1}}\bm{b}_{1}}_{\hat{\bm{a}}_{2}},\qquad$
$\displaystyle(\text{Combination view})$

where the first equation shows $\bm{b}_{2}$ is a multiplication of a matrix
and $\bm{a}_{2}$, i.e., project $\bm{a}_{2}$ onto the orthogonal complement
space of $\mathcal{C}{([\bm{b}_{1}])}$. The second equation shows $\bm{a}_{2}$
is a combination of $\bm{b}_{1}$ and $\bm{b}_{2}$. Clearly, the space spanned
by $\bm{b}_{1},\bm{b}_{2}$ is the same space spanned by
$\bm{a}_{1},\bm{a}_{2}$. The siguation is shown in Figure 11.2(a) in which we
choose the direction of $\bm{b}_{1}$ as the $x$-axis in the Cartesian
coordinate system. $\hat{\bm{a}}_{2}$ is the projection of $\bm{a}_{2}$ onto
the line $\bm{b}_{1}$. It can be clearly shown that the part of $\bm{a}_{2}$
perpendicular to $\bm{b}_{1}$ is $\bm{b}_{2}=\bm{a}_{2}-\hat{\bm{a}}_{2}$ from
the figure.

For the third vector $\bm{b}_{3}$, it must be perpendicular to both the
$\bm{b}_{1}$ and $\bm{b}_{2}$ which is actually the vector $\bm{a}_{3}$
subtracting its projection along the plane spanned by $\bm{b}_{1}$ and
$\bm{b}_{2}$

$\displaystyle\bm{b}_{3}$
$\displaystyle=\bm{a}_{3}-\frac{\bm{b}_{1}\bm{b}_{1}^{\top}}{\bm{b}_{1}^{\top}\bm{b}_{1}}\bm{a}_{3}-\frac{\bm{b}_{2}\bm{b}_{2}^{\top}}{\bm{b}_{2}^{\top}\bm{b}_{2}}\bm{a}_{3}=(\bm{I}-\frac{\bm{b}_{1}\bm{b}_{1}^{\top}}{\bm{b}_{1}^{\top}\bm{b}_{1}}-\frac{\bm{b}_{2}\bm{b}_{2}^{\top}}{\bm{b}_{2}^{\top}\bm{b}_{2}})\bm{a}_{3}\qquad$
$\displaystyle(\text{Projection view})$ (11.1)
$\displaystyle=\bm{a}_{3}-\underbrace{\frac{\bm{b}_{1}^{\top}\bm{a}_{3}}{\bm{b}_{1}^{\top}\bm{b}_{1}}\bm{b}_{1}}_{\hat{\bm{a}}_{3}}-\underbrace{\frac{\bm{b}_{2}^{\top}\bm{a}_{3}}{\bm{b}_{2}^{\top}\bm{b}_{2}}\bm{b}_{2}}_{\bar{\bm{a}}_{3}},\qquad$
$\displaystyle(\text{Combination view})$

where the first equation shows $\bm{b}_{3}$ is a multiplication of a matrix
and $\bm{a}_{3}$, i.e., project $\bm{a}_{3}$ onto the orthogonal complement
space of $\mathcal{C}{([\bm{b}_{1},\bm{b}_{2}])}$. The second equation shows
$\bm{a}_{3}$ is a combination of $\bm{b}_{1},\bm{b}_{2},\bm{b}_{3}$. We will
see this property is essential in the idea of QR decomposition. Again, it can
be shown that the space spanned by $\bm{b}_{1},\bm{b}_{2},\bm{b}_{3}$ is the
same space spanned by $\bm{a}_{1},\bm{a}_{2},\bm{a}_{3}$. The situation is
shown in Figure 11.2(b), in which we choose the direction of $\bm{b}_{2}$ as
the $y$-axis of the Cartesian coordinate system. $\hat{\bm{a}}_{3}$ is the
projection of $\bm{a}_{3}$ onto the line $\bm{b}_{1}$, $\bar{\bm{a}}_{3}$ is
the projection of $\bm{a}_{3}$ onto the line $\bm{b}_{2}$. It can be shown
that the part of $\bm{a}_{3}$ perpendicular to both $\bm{b}_{1}$ and
$\bm{b}_{2}$ is $\bm{b}_{3}=\bm{a}_{3}-\hat{\bm{a}}_{3}-\bar{\bm{a}}_{3}$ from
the figure.

Finally, we normalize each vector by dividing their length which produces
three orthonormal vectors $\bm{q}_{1}=\frac{\bm{b}_{1}}{||\bm{b}_{1}||}$,
$\bm{q}_{2}=\frac{\bm{b}_{2}}{||\bm{b}_{2}||}$,
$\bm{q}_{2}=\frac{\bm{b}_{2}}{||\bm{b}_{2}||}$.

(a) Project $\bm{a}_{2}$ onto the space perpendicular to $\bm{b}_{1}$.

(b) Project $\bm{a}_{3}$ onto the space perpendicular to
$\bm{b}_{1},\bm{b}_{2}$.

Figure 11.2: Gram-Schmidt process.

This idea can be extended to a set of vectors rather than only three. And we
call this process as Gram-Schmidt process. After this process, matrix $\bm{X}$
will be triangularized. The method is named after Jørgen Pedersen Gram and
Erhard Schmidt, but it appeared earlier in the work of Pierre-Simon Laplace in
the theory of Lie group decomposition.

The Gram–Schmidt process is not the only algorithm for QR decomposition.
Several other QR decomposition algorithms exist such as Householder
reflections and Givens rotations which are more reliable in the presence of
round-off errors. These QR decomposition methods may also change the order in
which the columns of $\bm{X}$ are processed.

####  B.4 Orthogonal vs Orthonormal

The vectors $\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{p}\in^{n}$ are orthogonal
when their dot products $\bm{q}_{i}^{\top}\bm{q}_{j}$ are zero whenever $i\neq
j$. When each vector is divided by its length, the vectors become orthogonal
unit vectors. Then the vectors $\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{p}$ are
orthonormal. We put the orthonormal vectors into a matrix $\bm{Q}$.

When $n\neq p$: the matrix $\bm{Q}$ is easy to work with because
$\bm{Q}^{\top}\bm{Q}=\bm{I}\in^{p\times p}$.

When $n=p$: the matrix $\bm{Q}$ is square, $\bm{Q}^{\top}\bm{Q}=\bm{I}$ means
that $\bm{Q}^{\top}=\bm{Q}^{-1}$, i.e., the transpose of $\bm{Q}$ is the
inverse of $\bm{Q}$. Then we also have $\bm{Q}\bm{Q}^{\top}=\bm{I}$, i.e.,
$\bm{Q}^{\top}$ is the two-sided inverse of $\bm{Q}$. We call this $\bm{Q}$ an
orthogonal matrix. Note here we use the term orthogonal matrix to mean the
matrix $\bm{Q}$ has orthonormal columns. The term orthonormal matrix is not
used for historical reasons.

####  B.5 Property of the QR Decomposition

Given any matrix $\bm{X}$, we have the property: $\mathcal{N}(\bm{X}^{\top})$
is the orthogonal complement of the column space $\mathcal{C}(\bm{X})$ in n:
$dim(\mathcal{N}(\bm{X}^{\top}))+dim(\mathcal{C}(\bm{X}))=n$; This is called
rank-nullity theorem. And the proof can be found in Appendix B (p. B). In
specific, from QR decomposition, we can find a basis for this subspace. In
singular value decomposition, we will also find the basis for
$\mathcal{N}(\bm{X})$ and $\mathcal{C}(\bm{X}^{\top})$.

######  Lemma 11.6 (Orthonormal Basis in n)

For full QR decomposition, we have the following property:

$\bullet$ $\\{\bm{q}_{1},\bm{q}_{2}\ldots,\bm{q}_{p}\\}$ is an orthonormal
basis of $\mathcal{C}(\bm{X})$;

$\bullet$ $\\{\bm{q}_{p+1},\bm{q}_{p+2},\ldots,\bm{q}_{n}\\}$ is an
orthonormal basis of $\mathcal{N}(\bm{X}^{\top})$.

Proof [of Lemma 11.6] From the Gram-Schmidt process, it is trivial that
$span\\{\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{k}\\}$ is equal to
$span\\{\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{k}\\}$ for all
$k\in\\{1,2,\ldots,p\\}$. Thus
$\mathcal{C}(\bm{X})=span\\{\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{p}\\}=span\\{\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{p}\\}$,
and $\\{\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{p}\\}$ is an orthonormal basis
for the column space of $\bm{X}$. As
$\mathcal{N}(\bm{X}^{\top})\bot\mathcal{C}(\bm{X})$,
$dim(\mathcal{N}(\bm{X}^{\top}))=n-dim(\mathcal{C}(\bm{X}))=n-p$. And the
space spanned by $\\{\bm{q}_{p+1},\bm{q}_{p+2},\ldots,\bm{q}_{n}\\}$ is also
$\bot\mathcal{C}(\bm{X})$ with dimension $n-p$. Thus,
$\\{\bm{q}_{p+1},\bm{q}_{p+2},\ldots,\bm{q}_{n}\\}$ is an orthonormal basis
for $\mathcal{N}(\bm{X}^{\top})$.

####  B.6 Computing the Reduced QR Decomposition via the Gram-Schmidt Process

We write out this form of QR Decomposition

$\bm{X}=\left[\begin{matrix}\bm{x}_{1}&\bm{x}_{2}&\ldots&\bm{x}_{p}\end{matrix}\right]=\left[\begin{matrix}\bm{q}_{1}&\bm{q}_{2}&\ldots&\bm{q}_{p}\end{matrix}\right]\begin{bmatrix}r_{11}&r_{12}&\dots&r_{1p}\\\
&r_{22}&\dots&r_{2p}\\\ &&\ddots&\vdots\\\
<EMAIL_ADDRESS><EMAIL_ADDRESS>
The orthogonal matrix $\bm{Q}$ can be easily calculated by the Gram-Schmidt
process. To see why we have the upper triangular matrix $\bm{R}$, we write out
these equations:

$\displaystyle\bm{x}_{1}$ $\displaystyle=r_{11}\bm{q}_{1}$
$\displaystyle=\sum_{i=1}^{1}r_{i1}\bm{q}_{1},$ $\displaystyle\bm{x}_{2}$
$\displaystyle=r_{12}\bm{q}_{1}+r_{22}\bm{q}_{2}$
$\displaystyle=\sum_{i=1}^{2}r_{i2}\bm{q}_{2},$ $\displaystyle\vdots$
$\displaystyle\bm{x}_{p}$
$\displaystyle=r_{1p}\bm{q}_{1}+r_{2p}\bm{q}_{2}+\ldots+r_{pp}\bm{q}_{p}$
$\displaystyle=\sum_{i=1}^{p}r_{ip}\bm{q}_{p},$

which coincide with the second equation of Equation (11.1) and make the shape
of an upper triangular matrix $\bm{R}$. And if we extend the idea of Equation
(11.1) into the $p$-th term, we get

$\displaystyle\bm{x}_{p}$
$\displaystyle=\sum_{i=1}^{p-1}(\bm{q}_{i}^{\top}\bm{x}_{p})\bm{q}_{i}+\bm{b}_{p}$
$\displaystyle=\sum_{i=1}^{p-1}(\bm{q}_{i}^{\top}\bm{x}_{p})\bm{q}_{i}+||\bm{b}_{p}||\bm{q}_{p},$

which implies we can gradually orthonormalize $\bm{X}$ to an orthonormal set
$\bm{Q}=[\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{p}]$ as in Algorithm 2.

Algorithm 2 Reduced QR Decomposition

1:matrix $\bm{X}$ has independent columns with size $n\times p$ and $n\geq p$;
$\displaystyle 1:\bm{q}_{1}$ $\displaystyle=\bm{x}_{1}/r_{11},$ $\displaystyle
r_{11}$ $\displaystyle=||\bm{x}_{1}||,$ $\displaystyle 2:\bm{q}_{2}$
$\displaystyle=(\bm{x}_{2}-r_{12}\bm{q}_{1})/r_{22},$ $\displaystyle r_{12}$
$\displaystyle=\bm{q}_{1}^{\top}\bm{x}_{2},r_{22}=||\bm{x}_{2}-r_{12}\bm{q}_{1}||,$
$\displaystyle 3:\bm{q}_{3}$
$\displaystyle=(\bm{x}_{3}-r_{13}\bm{q}_{1}-r_{12}\bm{q}_{2})/r_{33},$
$\displaystyle r_{13}$
$\displaystyle=\bm{q}_{1}^{\top}\bm{x}_{3},r_{23}=\bm{q}_{2}^{\top}\bm{x}_{3},r_{33}=||\bm{x}_{3}-r_{13}\bm{q}_{1}-r_{23}\bm{q}_{2}||,$
$\displaystyle\vdots$ $\displaystyle=\vdots$ $\displaystyle\vdots$
$\displaystyle n:\bm{q}_{p}$
$\displaystyle=(\bm{x}_{p}-\sum_{i=1}^{p-1}r_{ip}\bm{q}_{i})/r_{pp},$
$\displaystyle r_{ip}$
$\displaystyle=\bm{q}_{i}^{\top}\bm{x}_{p},r_{pp}=||\bm{x}_{p}-\sum_{i=1}^{p-1}r_{ip}\bm{q}_{i}||,\forall
i\in\\{1,2,\ldots,p-1\\},$

It can be shown Algorithm 2 requires $\sim 2np^{2}$ flops to compute a reduced
QR decomposition of an $n\times p$ matrix with independent columns and $n\geq
p$. See Lu (2021a) for more details on the complexity of QR decomposition.

####  B.7 Full QR Decomposition via the Gram-Schmidt Process

A full QR decomposition of an $n\times p$ matrix with independent columns goes
further by appending additional $n-p$ orthonormal columns to $\bm{Q}$ so that
it becomes an $n\times n$ orthogonal matrix. In addition, rows of zeros are
appended to $\bm{R}$ so that it becomes an $n\times p$ upper triangular
matrix. We refer to the additional columns in $\bm{Q}$ as silent columns and
additional rows in $\bm{R}$ as silent rows. The comparison of reduced QR
decomposition and full QR decomposition is shown in Figure 1.4 (p. 1.4) where
silent columns in $\bm{Q}$ are denoted in grey, blank entries are zero and
blue entries are elements that are not necessarily zero.

####  B.8 Dependent Columns

Previously, we assumed matrix $\bm{X}$ has independent columns. However, this
is not always necessary. Suppose in step $k$ of Algorithm 2, $\bm{x}_{k}$ is
in the plane spanned by $\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{k-1}$ which is
equivalent to the space spanned by
$\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{k-1}$, i.e., vectors
$\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{k}$ are dependent. Then $r_{kk}$ will be
zero and $\bm{q}_{k}$ does not exist because of the zero division. At this
moment, we simply pick $\bm{q}_{k}$ arbitrarily to be any normalized vector
that is orthogonal to
$\mathcal{C}([\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{k-1}])$ and continue the
Gram-Schmidt process. Again, for matrix $\bm{X}$ with dependent columns, we
have both reduced and full QR decomposition algorithms. We reformulate the
$k$-th step in the algorithm as follows:

$\bm{q}_{k}=\left\\{\begin{aligned}
&(\bm{x}_{k}-\sum_{i=1}^{k-1}r_{ik}\bm{q}_{i})/r_{kk},\qquad
r_{ik}=\bm{q}_{i}^{\top}\bm{x}_{k},r_{kk}=||\bm{x}_{k}-\sum_{i=1}^{k-1}r_{ik}\bm{q}_{i}||,&\mathrm{if\,}r_{kk}\neq
0,\\\
&\mathrm{pick\,one\,in\,}\mathcal{C}^{\bot}([\bm{q}_{1},\bm{q}_{2},\ldots,\bm{q}_{k-1}]),\qquad&\mathrm{if\,}r_{kk}=0.\end{aligned}\right.$

This idea can be further extended that when $\bm{q}_{k}$ does not exist, we
just skip the current steps. And add the silent columns in the end. In this
sense, QR decomposition for a matrix with dependent columns is not unique.
However, as long as you follow a systematic process or methodical procedure,
QR decomposition for any matrix is unique.

This finding can also help to decide whether a set of vectors are independent
or not. Whenever $r_{kk}$ in Algorithm 2 is zero, we report the vectors
$\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{k}$ are dependent and stop the algorithm
for “independent checking”.

####  B.9 Existence of the QR Decomposition via the Householder Reflector

In this section, we provide another view to find the QR decomposition via the
Householder reflector. We first give the formal definition of a Householder
reflector and we will discuss its properties.

######  Definition 11.1 (Householder Reflector)

Let $\bm{u}\in^{p}$ be a vector of unit length (i.e., $||\bm{u}||=1$). Then
$\bm{H}=\bm{I}-2\bm{u}\bm{u}^{\top}$ is said to be a Householder reflector,
a.k.a., Householder transformation. We call this $\bm{H}$ the Householder
reflector associated with the unit vector $\bm{u}$.

Then we have the following corollary from this definition.

######  Corollary 11.7 (Unreflected by Householder)

Any vector $\bm{v}$ that is perpendicular to $\bm{u}$ is left unchanged by the
Householder transformation.

The proof is trivial that
$(\bm{I}-2\bm{u}\bm{u}^{\top})\bm{v}=\bm{v}-2\bm{u}\bm{u}^{\top}\bm{v}=\bm{v}$.

Suppose $\bm{u}$ is a unit vector, and a vector $\bm{v}$ is perpendicular to
$\bm{u}$. Then any vector $\bm{x}$ on the plane can be decomposed into two
parts $\bm{x}=\bm{x}_{\bm{v}}+\bm{x}_{\bm{u}}$: the first one
$\bm{x}_{\bm{u}}$ is parallel to $\bm{u}$ and the second one $\bm{x}_{\bm{v}}$
is perpendicular to $\bm{u}$ (i.e., parallel to $\bm{v}$). Following from
Section B.1 (p. B.1), $\bm{x}_{\bm{u}}$ can be computed by
$\bm{x}_{\bm{u}}=\frac{\bm{u}\bm{u}^{\top}}{\bm{u}^{\top}\bm{u}}\bm{x}=\bm{u}\bm{u}^{\top}\bm{x}$.
We then transform this $\bm{x}$ by the Householder reflector associated with
$\bm{u}$,
$\bm{H}\bm{x}=(\bm{I}-2\bm{u}\bm{u}^{\top})(\bm{x}_{\bm{v}}+\bm{x}_{\bm{u}})=\bm{x}_{\bm{v}}-\bm{u}\bm{u}^{\top}\bm{x}=\bm{x}_{\bm{v}}-\bm{x}_{\bm{u}}$.
That is, the space perpendicular to $\bm{u}$ acts as a mirror and any vector
$\bm{x}$ is reflected by the Householder reflector associated with $\bm{u}$.
The situation is shown in Figure 11.3.

Figure 11.3: Demonstration of Householder reflector

If we know two vectors are reflected to each other, the next corollary tells
us how to find the associated Householder reflector.

######  Corollary 11.8 (Finding the Householder Reflector)

Suppose $\bm{x}$ is reflected to $\bm{y}$ by a Householder reflector, then the
Householder reflector is

$\bm{H}=\bm{I}-2\bm{u}\bm{u}^{\top},\text{ where
}\bm{u}=\frac{\bm{x}-\bm{y}}{||\bm{x}-\bm{y}||}.$

###### Remark 11.9

If $\bm{H}$ is a Householder reflector, then it has the following properties:

$\bullet$ $\bm{H}\bm{H}=\bm{I}$;

$\bullet$ $\bm{H}=\bm{H}^{\top}$;

$\bullet$ $\bm{H}^{\top}\bm{H}=\bm{H}\bm{H}^{\top}=\bm{I}$ such that
Householder reflector is an orthogonal matrix.

We note from the Gram-Schmidt section that QR decomposition is to use an
orthogonal to triangularize a matrix $\bm{X}$. The further idea is that, if we
have a set of orthogonal matrices that can make $\bm{X}$ to be triangular step
by step, then we can also recover the QR decomposition. Specifically, if we
have an orthogonal matrix $\bm{Q}_{1}$ that can introduce zeros to the 1-st
column of $\bm{X}$ except the entry (1,1), and an orthogonal matrix
$\bm{Q}_{2}$ that can introduce zeros to the 2-nd column except the entries
(2,1), (2,2), $\ldots$. Then, we can also find the QR decomposition.

For the way to introduce zeros, we could reflect the columns of the matrix to
a basis vector $\bm{e}_{1}$ whose entries are all zero except the first entry.