text:
# Iterative solvers and preconditioning

A successful family of methods, usually referred to as Conjugate
Gradient-like algorithms, or Krylov subspace methods, can be viewed as
Galerkin or least-squares methods applied to a linear system $Ax =b$.
This view is different from the standard approaches to deriving the
classical Conjugate Gradient method in the literature. Nevertheless,
the fundamental ideas of least squares and Galerkin approximations
can be used to derive the most
popular and successful methods for linear systems, and this is the
topic of the present chapter.  Such a view may increase the general
understanding of variational methods and their applicability.

Given a linear system

$$
\begin{equation}
Ax = b,\quad x,b\in\mathbb{R}^n,\ A\in\mathbb{R}^{n,n}
\label{_auto1} \tag{1}
\end{equation}
$$

and a start vector $x^0$,
we want to construct an iterative solution method that produces approximations
$x^1,x^2,\ldots$, which hopefully converge to the exact solution $x$.
In iteration no. $k$ we seek an approximation

$$
\begin{equation}
x^{k+1} = x^{k} + u,\quad u = \sum_{j=0}^k c_j{q}_j
\label{linalg:xk:update} \tag{2}
\end{equation}
$$

where ${q}_j\in\mathbb{R}^n$ are known vectors

and $c_j$ are constants to be determined.
To be specific, let $q_0,\ldots,q_k$ be basis vectors for $V_{k+1}$:

$$
V_{k+1} = \hbox{span}\{q_0,\ldots,q_k\}{\thinspace .}
$$

The associated inner product $(\cdot,\cdot )$ is here
the standard Euclidean inner product on $\mathbb{R}^n$.

The corresponding error in the equation $Ax =b$, the residual, becomes

$$
r^{k+1} = b -Ax^{k+1} =
r^{k} -\sum_{j=0}^kc_jA{q}_j {\thinspace .}
$$

# Conjugate gradient-like iterative methods

## The Galerkin method

Galerkin's method states that the error in the equation, the residual,
is orthogonal to the space $V_{k+1}$ where we seek the approximation to
the problem.

The Galerkin (or projection) method aims at finding $u=\sum_jc_jq_j\in V_{k+1}$
such that

$$
(r^{k+1},v)=0,\quad\forall v\in V_{k+1}{\thinspace .}
$$

This statement is
equivalent to the residual being orthogonal to each basis vector:

$$
\begin{equation}
(r^{k+1},{q}_i )=0,\quad i=0,\ldots,k{\thinspace .}
\label{linalg2:cg:eq1} \tag{3}
\end{equation}
$$

Inserting the expression for $r^{k+1}$ in ([3](#linalg2:cg:eq1))
gives a linear system
for $c_j$:

$$
\begin{equation}
\sum_{j=0}^k (A{q}_i,{q}_j )c_j = (r^{k},{q}_i),\quad i=0,\ldots,k
{\thinspace .}
\label{_auto2} \tag{4}
\end{equation}
$$

## The least squares method

The idea of the least-squares method is to
minimize the square of the norm of the
residual with respect to the free parameters $c_0,\ldots,c_k$.
That is, we minimize $(r^{k+1},r^{k+1})$:

$$
{\partial\over\partial c_i} (r^{k+1},r^{k+1}) =
2({\partial r^{k+1}\over
\partial c_i},r^{k+1}) =0,\quad i=0,\ldots,k{\thinspace .}
$$

Since $\partial r^{k+1} /\partial c_i = -A{q}_i$, this approach leads to
the following linear system:

$$
\begin{equation}
\sum_{j=0}^k (A{q}_i,A{q}_j)c_j = (r^{k+1},A{q}_i),\quad i=0,\ldots,k{\thinspace .}
\label{ch:linalg:ls:eq1} \tag{5}
\end{equation}
$$

## Krylov subspaces

To obtain a complete algorithm, we need to establish a rule to update
the basis $\mathcal{B}=\{q_0,\ldots,{q}_k\}$ for the next iteration.
That is, we need to
compute a new basis vector ${q}_{k+1}\in V_{k+2}$ such that

$$
\begin{equation}
\mathcal{B} =\{ {q}_0,\ldots,{q}_{k+1}\}
\label{linalg2:cg:Basis} \tag{6}
\end{equation}
$$

is a basis for the space $V_{k+2}$ that is used in the next iteration.
The present family of methods applies the *Krylov space*, where $V_k$
is defined as

$$
\begin{equation}
V_{k} = \mbox{span} \left\lbrace r^0,Ar^0,A^2r^0,\ldots
A^{k-1}r^0 \right\rbrace {\thinspace .}
\label{_auto3} \tag{7}
\end{equation}
$$

Some frequent names
of the associated iterative methods are therefore
{Krylov subspace iterations}, Krylov projection methods,
or simply Krylov methods.
It is a fact that $V_{k} \subset V_{k+1}$ and that $r^0
Ar^0,\ldots A^{k-1}r^0$ are linearly independent
vectors.

## Computation of the basis vectors

A potential formula for updating ${q}_{k+1}$, such that
${q}_{k+1}\in V_{k+2}$, is

$$
\begin{equation}
{q}_{k+1} =  r^{k+1} + \sum_{j=0}^k\beta_j{q}_j{\thinspace .}
\label{linalg:q:update1} \tag{8}
\end{equation}
$$

(Since $r^{k+1}$ involves $Aq_k$, and $q_k\in V_{k+1}$, multiplying
by $A$ raises the dimension of the Krylov space by 1, so $Aq_k\in V_{k+2}$.)
The free parameters $\beta_j$ can be used to enforce desirable
orthogonality properties of ${q}_0,\ldots,{q}_{k+1}$. For example,
it is convenient to require that the coefficient matrices in the linear
systems for $c_0,\ldots,c_k$ are diagonal.
Otherwise, we must solve a $(k+1)\times (k+1)$ linear system in each iteration.
If $k$ should approach
$n$, the systems for the coefficients $c_i$ are of
the same size as our original system $Ax =b$!
A diagonal matrix, however, ensures an efficient closed form solution for
$c_0,\ldots,c_k$.

To obtain a diagonal coefficient matrix, we require in Galerkin's method
that

$$
(A{q}_i,{q}_j)=0\quad \mbox{when}\ i\neq j,
$$

whereas we in the least-squares method require

$$
(A{q}_i,A{q}_j)=0\quad \mbox{when}\ i\neq j{\thinspace .}
$$

We can define
the inner product

$$
\begin{equation}
\langle u,v\rangle\equiv (Au,v ) =u^TAv,
\label{_auto4} \tag{9}
\end{equation}
$$

provided $A$ is symmetric and positive definite. Another
useful inner product is

$$
\begin{equation}
[u,v ]\equiv (Au,Av) = u^TA^TAv {\thinspace .}
\label{_auto5} \tag{10}
\end{equation}
$$

These inner products will be referred to as the $A$ product, with the
associated $A$ norm, and the $A^TA$ product, with the associated
$A^TA$ norm.

The orthogonality condition on the ${q}_i$ vectors are then
$\langle {q}_{k+1},{q}_i\rangle =0$ in the Galerkin method
and  $[{q}_{k+1},{q}_i]=0$ in the least-squares method, where $i$
runs from 0 to $k$.
A standard Gram-Schmidt process can be used for constructing
${q}_{k+1}$ orthogonal to ${q}_0,\ldots,{q}_k$. This leads to the determination
of the $\beta_0,\ldots,\beta_k$ constants as

$$
\begin{equation}
\beta_i = { \langle r^{k+1},{q}_i\rangle\over\langle{q}_i,{q}_i\rangle}
\quad\hbox{(Galerkin)} \label{linsys:betaiG} \tag{11}
\end{equation}
$$

$$
\begin{equation}
\beta_i = { [r^{k+1},{q}_i]\over [{q}_i,{q}_i] }
\quad\hbox{(least squares)}
\label{_auto6} \tag{12}
\end{equation}
$$

for $i=0,\ldots,k$.

## Computation of a new solution vector

The orthogonality condition on the basis vectors ${q}_i$ leads to
the following solution for $c_0,\ldots,c_k$:

$$
\begin{equation}
c_i = { (r^{k},{q}_i)\over \langle {q}_i,{q}_i\rangle}
\quad\hbox{(Galerkin)} \label{linsys:alpha:G} \tag{13}
\end{equation}
$$

$$
\begin{equation}
c_i = { (r^{k},A{q}_i)\over [ {q}_i,{q}_i]}
\quad\hbox{(least squares)} \label{linsys:alpha:LS} \tag{14}
\end{equation}
$$

In iteration $k$, $(r^{k},{q}_i)=0$ and
$(r^{k},A{q}_i)=0$, for $i=0,\ldots,k-1$, in the Galerkin and
least squares case, respectively.  Hence, $c_i =0$, for $i=0,\ldots,
k-1$. In other words,

$$
x^{k+1} = x^{k} + c_k{q}_k {\thinspace .}
$$

When $A$ is symmetric and positive definite, one can show that also
$\beta_i=0$, for $0=1,\ldots,k-1$, in both the Galerkin and least
squares methods.  This means that $x^k$ and
${q}_{k+1}$ can be updated using only ${q}_k$ and not the previous
${q}_0,\ldots,{q}_{k-1}$ vectors.  This property has of course dramatic
effects on the storage requirements of the algorithms as the number of
iterations increases.

For the suggested algorithms to work, we must require that the
denominators in ([13](#linsys:alpha:G)) and ([14](#linsys:alpha:LS)) do
not vanish.  This is always fulfilled for the least-squares method,
while a (positive or negative) definite matrix $A$ avoids break-down
of the Galerkin-based iteration (provided ${q}_i \neq 0$).

The Galerkin solution method for linear systems was originally devised
as a *direct method* in the 1950s.  After $n$ iterations the exact
solution is found in exact arithmetic, but at a higher cost than when
using Gaussian elimination. Naturally, the method did not receive
significant popularity before researchers discovered (in the beginning
of the 1970s) that the method could produce a good approximation to
$x$ for $k\ll n$ iterations.  The method, called the Conjugate
Gradient method, has from then on caught considerable interest as an
iterative scheme for solving linear systems arising from PDEs
discretized such that the coefficient matrix becomes sparse.

Finally, we mention how to terminate the iteration.
The simplest criterion is $||r^{k+1}||\leq\epsilon_r$, where
$\epsilon_r$ is a small prescribed tolerance.
Sometimes it is more appropriate to use a relative residual,
$||r^{k+1}||/||r^0||\leq\epsilon_r$.
Termination criteria for Conjugate Gradient-like methods is a subject
on its own.

## Summary of the least squares method

In the algorithm below, we have summarized
the computational steps in the least-squares method.
Notice that we update the residual recursively instead of using
$r^k=b - Ax^k$ in each iteration since we then avoid a possibly
expensive matrix-vector product.

1. given a start vector $x^0$,
   compute $r^0=b - Ax^0$ and set ${q}_0 =r^0$.

2. for $k=0,1,2,\ldots$ until termination criteria are fulfilled:

a. $c_k = {(r^{k},A{q}_k)/ [{q}_k,{q}_k]}$

b. $x^{k+1} = x^{k} + c_k{q}_k$

c. $r^{k+1} = r^{k} - c_kA{q}_k$

d. if $A$ is symmetric then

1. $\beta_k = {[r^{k+1},{q}_k]/ [{q}_k,{q}_k]}$

a. ${q}_{k+1} = r^{k+1} - \beta_k{q}_k$

b. else

1. $\beta_j = {[r^{k+1},{q}_j]/ [{q}_j,{q}_j]},\quad j=0,\ldots,k$

2. ${q}_{k+1} = r^{k+1} - \sum_{j=0}^k\beta_j{q}_j$

### Remark

The algorithm above is just a summary of the steps in the derivation
of the least-squares method and should not be directly used for
practical computations without further developments.

## Truncation and restart

When $A$ is nonsymmetric, the storage requirements of ${q}_0,\ldots,{q}_k$
may be prohibitively large. It has become a standard trick to
either *truncate* or *restart* the algorithm.
In the latter case one restarts the algorithm every $K+1$-th step, i.e.,
one aborts the iteration and starts the algorithm again with $x^0=x^K$.
The other alternative is to truncate the sum $\sum_{j=0}^k\beta_j{q}_j$
and use only the last $K$ vectors:

$$
x^{k+1} = x^{k} + \sum_{j=k-K}^k \beta_j{q}_j{\thinspace .}
$$

Both the restarted and truncated version of the algorithm require
storage of only $K+1$ basis vectors ${q}_{k-K},\ldots,{q}_k$.
The basis vectors are also often called *search direction vectors*.
The truncated version of the least-squares algorithm above
is widely known as *Generalized Minimum Residuals*, shortened as GMRES,
or GMRES$(K)$ to explicitly indicate the
number of search direction vectors.
In the literature one encounters the name
*Generalized Conjugate Residual method*, abbreviated CGR, for
the restarted version of Orthomin. When $A$ is symmetric, the method
is known under the name *Conjugate Residuals*.

## Summary of the Galerkin method

In case of Galerkin's method, we assume that $A$ is symmetric and
positive definite. The resulting computational procedure
is the famous Conjugate Gradient
method.
Since $A$ must be symmetric, the recursive update of
${q}_{k+1}$ needs only one previous search direction vector ${q}_k$, that is,
$\beta_j=0$ for $j<k$.

1. Given a start vector $x^0$,
   compute $r^0=b - Ax^0$ and set ${q}_0 =r^0$.

2. for $k=1,2,\ldots$ until termination criteria are fulfilled:

a. $c_k = {(r^{k},{q}_k) / \langle {q}_k,{q}_k\rangle}$

b. $x^{k} = x^{k-1} + c_k{q}_k$

c. $r^{k} = r^{k-1} - c_kA{q}_k$

d. $\beta_k = {\langle r^{k+1},{q}_k\rangle / \langle{q}_k,{q}_k\rangle}$

e. ${q}_{k+1} = r^{k+1} - \beta_k{q}_k$

The previous remark that the listed algorithm is just a summary of the
steps in the solution procedure, and not an efficient algorithm
that should be implemented in its present form, must be repeated here.
In general, we recommend to rely on some high-quality linear algebra library
that offers an implementation of the Conjugate gradient method.

**The computational nature of Conjugate gradient-like methods.**

Looking at the Galerkin and least squares algorithms above,
one notice that the matrix $A$ is only used in matrix-vector
products. This means that it is sufficient
to store only the nonzero entries of $A$.
The rest of the algorithms consists of vector operations of the
type $y \leftarrow ax + y$,
the slightly more general variant $q\leftarrow ax +y$, as well as
inner products.

## A framework based on the error

Let us define the error $e^k = x - x^k$. Multiplying this equation by $A$
leads to the well-known relation between the error and the residual
for linear systems:

$$
\begin{equation}
Ae^k = r^k {\thinspace .}
\label{linalg:erroreq} \tag{15}
\end{equation}
$$

Using $r^k = Ae^k$ we can reformulate the Galerkin and
least-squares methods in terms of the error.  The Galerkin method can
then be written

$$
\begin{equation}
(r^k,{q}_i ) = (Ae^k,{q}_i) = \langle e^k,{q}_i\rangle = 0,\quad
i=0,\ldots,k{\thinspace .}
\label{_auto7} \tag{16}
\end{equation}
$$

For the least-squares method we obtain

$$
\begin{equation}
(r^k,A{q}_i) = [e^k,{q}_i] =0,\quad i=0,\ldots,k{\thinspace .}
\label{_auto8} \tag{17}
\end{equation}
$$

This means that

$$
\begin{align*}
\langle e^k, v\rangle  &= 0\quad\forall v\in V_{k+1}\hbox{ (Galerkin)}\\
\lbrack e^k, v \rbrack  &= 0\quad\forall v\in V_{k+1}\hbox{ (least-squares)}
\end{align*}
$$

In other words, the error is $A$-orthogonal to the space $V_{k+1}$ in the
Galerkin method, whereas the error is $A^TA$-orthogonal to $V_{k+1}$ in
the least-squares method. When the error is orthogonal to a space, we
find the best approximation in the associated norm to the solution in
that space. Specifically here, it means that for a symmetric and positive
definite $A$, the Conjugate gradient method finds the optimal adjustment
in $V_{k+1}$ of the vector $x^k$ (in the $A$-norm)
in the update for $x^{k+1}$. Similarly, the
least squares formulation finds the optimal adjustment in $V_{k+1}$
measured in the $A^TA$-norm.

A lot of Conjugate gradient-like methods were developed in the 1980s and
1990s, some of the most popular methods do not fit
directly into the framework presented here.  The theoretical
similarities between the methods are covered in the literature, and
we refer to it for algorithms and practical
comments related to widespread methods, such as the SYMMLQ method (for
symmetric indefinite systems), the Generalized Minimal Residual
(GMRES) method, the BiConjugate Gradient (BiCG) method, the
Quasi-Minimal Residual (QMR) method, and the BiConjugate Gradient
Stabilized (BiCGStab) method.  When $A$ is symmetric and positive
definite, the Conjugate gradient method is the optimal choice with
respect to computational efficiency, but when $A$ is nonsymmetric,
the performance of the methods is strongly problem dependent, and one
needs to experiment.

# Preconditioning

## Motivation and Basic Principles

The Conjugate Gradient method has been subject to extensive analysis,
and its convergence properties are well understood.
To reduce the initial error $e^0 =x -x^0$ with a factor
$0 <\epsilon\ll 1$ after $k$ iterations, or more precisely,
$||e^k||_{A}\leq\epsilon ||e^0||_{A}$, it can be shown that
$k$ is bounded by

$$
{1\over2}\ln{2\over\epsilon}\sqrt{\kappa},
$$

where $\kappa$ is the ratio of the largest and smallest eigenvalue of
$A$. The quantity $\kappa$  is commonly referred to as
the spectral *condition number*.
Common finite element and finite difference discretizations of
Poisson-like PDEs lead to $\kappa\sim h^{-2}$, where $h$ denotes the
mesh size. This implies that the Conjugate Gradient method converges
slowly in PDE problems with fine meshes, as the number of iterations is
proportional to $h^{-1}$.

To speed up the Conjugate Gradient method, we should
manipulate the eigenvalue distribution. For instance, we could reduce
the condition number $\kappa$. This can be achieved by so-called
*preconditioning*. Instead of applying the iterative method to
the system $Ax =b$, we multiply by a matrix $M^{-1}$ and
apply the iterative method to the mathematically equivalent system

$$
\begin{equation} M^{-1}Ax = M^{-1}b {\thinspace .} \label{_auto9} \tag{18}
\end{equation}
$$

The aim now is to construct a nonsingular *preconditioning matrix or algorithm*
such that $M^{-1}A$
has a more favorable condition number than $A$.
We remark that we use the notation $M^{-1}$ here to indicate that it should resemble the inverse
of the matrix $A$.

For increased flexibility we can write
$M^{-1} = C_LC_R$ and transform the system according to

$$
\begin{equation}
C_LAC_R y = C_Lb,\quad  y =C_R^{-1}x ,
\label{_auto10} \tag{19}
\end{equation}
$$

where $C_L$ is the *left* and $C_R$ is the *right*
preconditioner. If the original coefficient matrix $A$ is symmetric
and positive definite, $C_L = C_R^T$ leads to preservation of these
properties in the transformed system. This is important when applying
the Conjugate Gradient method to the preconditioned linear system.
Even if $A$ and $M$ are symmetric and positive definite,
$M^{-1}A$ does not necessarily inherit these properties.  It
appears that for practical purposes one can express the iterative
algorithms such that it is sufficient to work with a single
preconditioning matrix $M$ only.  We
shall therefore speak of preconditioning in terms of the left
preconditioner $M$ in the following.

## Use of the preconditioning matrix in the iterative methods

Optimal convergence for the Conjugate Gradient method is
achieved when the coefficient matrix $M^{-1}A$ equals the identity matrix
$I$ and only one iteration is required. In the algorithm we need to perform matrix-vector products
$M^{-1}Au$ for an arbitrary $u\in\mathbb{R}^n$.
This means that we have to solve a linear system with $M$ as coefficient
matrix in each iteration since we implement the product $y =M^{-1}Au$
in a two step fashion: First we compute $v = Au$ and then we
solve the linear system $My =v$ for $y$. The optimal choice $M =A$
therefore involves the solution of $Ay =v$ in each iteration, which
is a problem of the same complexity as our original system $Ax =b$.
The strategy must hence be to compute an $M^{-1} \approx A^{-1} $ such that the
algorithmic operations involved in the inversion of $M$ are cheap.

The preceding discussion motivates the following
demands on the preconditioning matrix $M$:

  * $M^{-1}$ should be a good approximation to $A$,

  * $M^{-1}$ should be inexpensive to compute,

  * $M^{-1}$ should be sparse in order to minimize storage requirements,

  * linear systems with $M$ as coefficient matrix must be efficiently solved.

Regarding the last property, such systems must be solved in
$\mathcal{O}(n)$ operations, that is, a complexity of the same order
as the vector updates in the Conjugate Gradient-like algorithms.
These four properties are contradictory and some sort of compromise
must be sought.

## Classical iterative methods as preconditioners

The simplest possible iterative method for solving $Ax=b$ is

$$
x^{k+1} = x^{k} + r^{k} {\thinspace .}
$$

Applying this method to the preconditioned system
$M^{-1}Ax =M^{-1}b$ results in the scheme

$$
x^{k+1} = x^{k} + M^{-1}r^{k} ,
$$

which is nothing but the formula for classical iterative methods such as
the Jacobi method, the Gauss-Seidel method, SOR (Successive over-relaxation),
and SSOR (Symmetric successive over-relaxation).
This motivates for choosing $M$ as one iteration of these classical methods.
In particular, these methods provide simple formulas for $M$:

  * Jacobi preconditioning: $M =D$.

  * Gauss-Seidel preconditioning: $M = D + L$.

  * SOR preconditioning: $M = \omega^{-1}D +L$.

  * SSOR preconditioning:  $M = (2-\omega )^{-1} \left( \omega^{-1}D + L\right) \left( \omega^{-1}D\right)^{-1} \left( \omega^{-1}D + U\right)$

Turning our attention to the four requirements of the preconditioning
matrix, we realize that the suggested $M$ matrices do not demand
additional storage, linear systems with $M$ as coefficient matrix are
solved effectively in $\mathcal{O}(n)$ operations, and $M$ needs no
initial computation.  The only questionable property is how well $M$
approximates $A$, and that is the weak point of using classical
iterative methods as preconditioners.

The Conjugate Gradient method can only utilize the Jacobi and SSOR
preconditioners among the classical iterative methods, because the $M$
matrix in that case is on the form $M^{-1}=C_LC_L^T$, which is
necessary to ensure that the coefficient matrix of the preconditioned
system is symmetric and positive definite.  For certain PDEs, like the
Poisson equation, it can be shown that the SSOR preconditioner reduces
the condition number with an order of magnitude, i.e., from
$\mathcal{O}(h^{-2})$ to $\mathcal{O}(h^{-1})$, provided we use the
optimal choice of the relaxation parameter $\omega$.  According to
experiments, however, the performance of the SSOR preconditioned
Conjugate Gradient method is not very sensitive to the choice of
$\omega$.

## Incomplete factorization preconditioners

Imagine that we choose $M =A$ and solve systems $My =v$ by a
direct method. Such methods typically first compute the LU
factorization $M =\bar L\bar U$ and thereafter perform two
triangular solves. The lower and upper triangular factors $\bar L$ and
$\bar U$ are computed from a Gaussian elimination procedure.
Unfortunately, $\bar L$ and $\bar U$ contain nonzero values, so-called
*fill-in*, in many locations where the original matrix $A$ contains
zeroes. This decreased sparsity of $\bar L$ and $\bar U$ increases
both the storage requirements and the computational efforts related to
solving systems $My =v$.  An idea to improve the situation is to
compute *sparse* versions of the factors $\bar L$ and $\bar U$. This
is achieved by performing Gaussian elimination, but neglecting the
fill-in (!). In this way, we can compute approximate factors $\widehat L$
and $\widehat U$ that become as sparse as $A$.  The storage
requirements are hence only doubled by introducing a preconditioner,
and the triangular solves become an $\mathcal{O}(n)$ operation since
the number of nonzeroes in the $\widehat L$ and $\widehat U$ matrices
(and $A$) is $\mathcal{O}(n)$ when the underlying PDE is discretized
by finite difference or finite element methods.  We call $M
=\widehat L\widehat U$ an *incomplete LU factorization*
preconditioner, often just referred to as the ILU preconditioner.

Instead of throwing away all fill-in entries, we can add them to the
main diagonal. This yields the *modified incomplete LU factorization*
(MILU).  One can also allow a certain amount of fill-in
to improve the quality of the preconditioner, at the expense of
more storage for the factors. Libraries offering ILU preconditioners
will then often have a tolerance for how small a fill-in element in
the Gaussian elimination must be to be dropped, and a parameter that
governs the maximum number of nonzeros per row in the factors.
A popular ILU implementation is the open source C code
[SuperLU](https://fs.hlrs.de/projects/craydoc/docs_merged/books/S-6532-30/S-6532-30.pdf). This preconditioner is available to Python programmers through the
`scipy.sparse.linalg.spilu` function.

The general algorithm for MILU preconditioning follows the steps of
traditional exact Gaussian elimination, except that
we restrict the computations to the nonzero entries in $A$.
The factors $\widehat L$ and $\widehat U$ can be stored directly in the
sparsity structure of $A$, that is, the algorithm overwrites a copy
$M$ of $A$ with
its MILU factorization. The steps in the MILU factorizations are
listed below.

1. Given a sparsity pattern as an index set $\mathcal{I}$,
   copy $M_{i,j}\leftarrow A_{i,j}$, $i,j=1,\ldots,n$

2. for $k=1,2,\ldots, n$

a. for $i=k+1,\ldots,n$

     * if $(i,k)\in\mathcal{I}$ then

a. $M_{i,k}\leftarrow M_{i,k}/M_{k,k}$

     * else

a. $M_{i,k} = 0$

     * $r=M_{i,k}$

     * for $j=k+1,\ldots,n$

       * if $j=i$ then

         * $M_{i,j}\leftarrow M_{i,j} - rM_{k,j} + \sum_{p=k+1}^n
           \left( M_{i,p} - rM_{k,p}\right)$

       * else

         * if $(i,j)\in\mathcal{I}$ then

           * $M_{i,j}\leftarrow M_{i,j} - rM_{k,j}$

         * else

           * $M_{i,j} =0$

We also remark here that the algorithm above needs careful refinement
before it should be implemented in a code. For example, one will not
run through a series of $(i,j)$ indices and test for each of them if
$(i,j)\in\mathcal{I}$. Instead one should run more directly through
the sparsity structure of $A$.

## Preconditioners developed for solving PDE problems

The above mentioned preconditioners are general purpose
preconditioners that may be applied to solve any linear system
and may speed up the solution process significantly. For linear
systems arising from PDE problems there are however often more
efficient preconditioners that exploit the fact that the matrices
come from discretized PDEs. These preconditioners often enable
the solution of linear systems involving millions of unknowns
in a matter of seconds or less on a regular desktop computer.
This means that solving these huge linear systems often
takes less time than printing them file. Moreover, these preconditioners
may be used on parallel computers in a scalable fashion such that
simulations involving billions of degrees of freedom are available
e.g. by cloud services or super-computing centers.

The most efficient preconditioners are the so-called multilevel preconditioners
(e.g. multigrid  and domain decomposition) which in particular for
Poisson type problems yields error reduction of a factor 10 per iteration
and hence a typical iteration count of the iterative method of 5-10.
These preconditioners utilize the fact that different parts of the
solution can be localized or represented on a coarser resolution during
the computations before being glued together to a close match to
the true solution. While a proper description of these methods
are beyond the scope here, we mention that black-box preconditioners
are implemented in open source frameworks such as PETSc, Hypre, and Trilinos
and these frameworks are used by most available finite element software (such as FEniCS).

While these multilevel preconditioners are available and efficient for PDEs similar
to the Poisson problem there is currently no black-box solver that handles
general systems of PDEs. We mention, however, a promising approach for
tackling some systems of PDEs, namely the so-called operator preconditioning
framework. This framework utilize tools from functional analysis to deduct
appropriate block preconditioners typically involving blocks composed
of for instance multilevel preconditioners in a systematic construction.
The framework has extended the use of Poisson type multilevel preconditioners
to many systems of PDEs with success.

max_stars_repo_path: 9- cg.ipynb

max_stars_repo_name: mbarzegary/finite-element-intro

max_issues_repo_path: 9- cg.ipynb

max_issues_repo_name: mbarzegary/finite-element-intro

max_forks_repo_path: 9- cg.ipynb

max_forks_repo_name: mbarzegary/finite-element-intro


text:
```python
import numpy as np
import sympy as sp
from sympy import Matrix, Eq
from sympy.physics.quantum import TensorProduct

a, b, c, d, e, f, g, h, i = sp.symbols('a b c d e f g h i')
A = Matrix([[a, b, c], [d, e, f], [g, h, i]])
A_T = A.transpose()
```

$A + A'$ is symmetric**

$(A + A')' \equiv A + A'$

```python
Eq((A + A_T).transpose(), A + A_T)
```

$\displaystyle \text{True}$

$A - A'$ is skew-symmetric

$(A - A')' \equiv -(A - A')$

```python
Eq((A - A_T).transpose(), A_T - A)
```

$\displaystyle \text{True}$

**Square matrix = Symmetric matrix + Anti-symmetric matrix**
$A = \frac{1}{2} (A + A') + \frac{1}{2} (A - A')$

```python
Eq(A, 1/2 * (A + A_T) + 1/2 * (A - A_T))
```

$\displaystyle \text{True}$

```python

```

```python
a = np.array([[9/2, 3], [15/2, 6]])

np.linalg.det(a)
```

    4.4999999999999964

```python
a, b, c, d, e, f, g, h, i = sp.symbols('a b c d e f g h i')
m1 = Matrix([[a,b],[c,d]])
m2 = Matrix([[e,f],[g,h]])
P = TensorProduct(m1, m2)
```

```python
m1
```

$\displaystyle \left[\begin{matrix}a & b\\c & d\end{matrix}\right]$

```python
m2
```

$\displaystyle \left[\begin{matrix}e & f\\g & h\end{matrix}\right]$

```python
P
```

$\displaystyle \left[\begin{matrix}a e & a f & b e & b f\\a g & a h & b g & b h\\c e & c f & d e & d f\\c g & c h & d g & d h\end{matrix}\right]$

```python
a, b, c, d, e, f, g, h, i = sp.symbols('a b c d e f g h i')
m1 = Matrix([[a],[b]])
m2 = Matrix([[c],[d]])
P = TensorProduct(m1, m2)
```

```python
m1
```

$\displaystyle \left[\begin{matrix}a\\b\end{matrix}\right]$

```python
m2
```

$\displaystyle \left[\begin{matrix}c\\d\end{matrix}\right]$

```python
P
```

$\displaystyle \left[\begin{matrix}a c\\a d\\b c\\b d\end{matrix}\right]$

```python
a, b, c, d, e, f, g, h, i = sp.symbols('a b c d e f g h i')
m1 = Matrix([a, b])
m2 = Matrix([c,d])
P = TensorProduct(m1, m2)
```

```python
m1
```

$\displaystyle \left[\begin{matrix}a\\b\end{matrix}\right]$

```python
m2
```

$\displaystyle \left[\begin{matrix}c\\d\end{matrix}\right]$

```python
P
```

$\displaystyle \left[\begin{matrix}a c\\a d\\b c\\b d\end{matrix}\right]$

```python

```

ext: ipynb

max_stars_repo_name: Valentine-Efagene/Jupyter-Notebooks

max_issues_repo_name: Valentine-Efagene/Jupyter-Notebooks

max_forks_repo_name: Valentine-Efagene/Jupyter-Notebooks


text:
# TIEGCM

The kamodofied tiegcm requires the pytiegcm reader available on github:

    pip install -e git+https://github.com/asher-pembroke/pytiegcm.git#egg=pytiegcm

The TIEGCM reader also requires `netcdf4` (and numpy, scipy, and pandas, but those are required by kamodo)

    pip install netcdf4

!!! warning
    This is a work in progress!

```python
filename = '/Users/apembrok/Work/TIEGCM/sample_data/jasoon_shim_040118_IT_1/s001.nc' # ccmc run (no zgmid)
fname = '/Users/apembrok/Downloads/2013.03.01.tie-gcm.data/s003.nc' # (no zgmid)
```

```python
from kamodo import readers
```

```python
from kamodo.readers.tiegcm import TIEGCM_Kamodo
```

By default, the kamodofied TIEGCM reader will load all 3d and 4d variables from the model. For performance reasons, you may want to only include certain variables for analysis. To do so, pass a list of variables in addition to the filename:

```python
tiegcm = TIEGCM_Kamodo(fname, ['TN', 'UN', 'VN', 'EFLUX', 'Z'])

```

    initializing tiegcm with /Users/apembrok/Downloads/2013.03.01.tie-gcm.data/s003.nc

    /Users/apembrok/git/kamodo/docs/notebooks/src/pytiegcm/tiegcm/tiegcm.py:213: UserWarning:

    WARNING: missing_value not used since it
    cannot be safely cast to variable data type

    registered 5 variables

```python
tiegcm
```

\begin{equation}\operatorname{TN}{\left(t,ilev,lat,lon \right)} [K] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{UN}{\left(t,ilev,lat,lon \right)} [cm/s] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{VN}{\left(t,ilev,lat,lon \right)} [cm/s] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}Z{\left(t,ilev,lat,lon \right)} [cm] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{EFLUX}{\left(t,lat,lon \right)} = \lambda{\left(t,lat,lon \right)}\end{equation}\begin{equation}\operatorname{ilev}{\left(t,z,lat,lon \right)} = \lambda{\left(t,z,lat,lon \right)}\end{equation}

Below is an attempt to use function composition to solve variables as a function of height.

```python
from scipy.interpolate import interp1d
```

```python
import numpy as np
z = np.squeeze(tiegcm.Z(tiegcm._time[0], lat = 0, lon = 0))
```

```python
tiegcm['T_N(t, z, lat, lon)'] = 'TN(t,ilev, lat, lon)'
```

```python
tiegcm
```

\begin{equation}\operatorname{TN}{\left(t,ilev,lat,lon \right)} [K] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{UN}{\left(t,ilev,lat,lon \right)} [cm/s] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{VN}{\left(t,ilev,lat,lon \right)} [cm/s] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}Z{\left(t,ilev,lat,lon \right)} [cm] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{EFLUX}{\left(t,lat,lon \right)} = \lambda{\left(t,lat,lon \right)}\end{equation}\begin{equation}\operatorname{ilev}{\left(t,z,lat,lon \right)} = \lambda{\left(t,z,lat,lon \right)}\end{equation}\begin{equation}\operatorname{T_{N}}{\left(t,z,lat,lon \right)} = \operatorname{TN}{\left(t,\operatorname{ilev}{\left(t,z,lat,lon \right)},lat,lon \right)}\end{equation}

```python
import numpy as np

@np.vectorize
def lev(self, t, z, lat, lon):
    # 1) meshgrid and squeeze the shapes of time, lat, lon if they are not the same
    # 2) calculate z_levels for all points
    # 3) construct rgi for (time, z_level, lat, lon)
    # 4) interpolate over z

    z_levels = self.Z(time = t, lat = lat, lon = lon)
    level = scipy.interpolate.interp1d(z_levels, ilev, bounds_error=False)
    return level(z)
tiegcm.lev = lev
```

```python
tiegcm['ilev'] = lev
```

```python
tiegcm['T_N(t, z, lat, lon)'] = 'TN(t, ilev, lat, lon)'
```

```python
tiegcm
```

\begin{equation}\operatorname{TN}{\left(t,ilev,lat,lon \right)} [K] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{UN}{\left(t,ilev,lat,lon \right)} [cm/s] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{VN}{\left(t,ilev,lat,lon \right)} [cm/s] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}Z{\left(t,ilev,lat,lon \right)} [cm] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{EFLUX}{\left(t,lat,lon \right)} = \lambda{\left(t,lat,lon \right)}\end{equation}\begin{equation}\operatorname{ilev}{\left(t,z,lat,lon \right)} = \lambda{\left(t,z,lat,lon \right)}\end{equation}\begin{equation}\operatorname{T_{N}}{\left(t,z,lat,lon \right)} = \operatorname{TN}{\left(t,\operatorname{ilev}{\left(t,z,lat,lon \right)},lat,lon \right)}\end{equation}

```python
tiegcm['eflux'] = 'EFLUX(t, lat, lon)'
```

```python
tiegcm
```

\begin{equation}\operatorname{TN}{\left(t,ilev,lat,lon \right)} [K] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{UN}{\left(t,ilev,lat,lon \right)} [cm/s] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{VN}{\left(t,ilev,lat,lon \right)} [cm/s] = \lambda{\left(t,ilev,lat,lon \right)}\end{equation}\begin{equation}\operatorname{EFLUX}{\left(t,lat,lon \right)} = \lambda{\left(t,lat,lon \right)}\end{equation}\begin{equation}t{\left(ilev \right)} = \lambda{\left(ilev \right)}\end{equation}\begin{equation}\operatorname{eflux}{\left(ilev,lat,lon \right)} = \operatorname{EFLUX}{\left(t{\left(ilev \right)},lat,lon \right)}\end{equation}

```python
%load_ext autoreload
%autoreload 2
```

    The autoreload extension is already loaded. To reload it, use:
      %reload_ext autoreload

```python

```

max_stars_repo_name: EnsembleGovServices/Kamodo-ccmc-readers

max_issues_repo_name: EnsembleGovServices/Kamodo-ccmc-readers

max_forks_repo_name: EnsembleGovServices/Kamodo-ccmc-readers


text:
$$\newcommand{\ket}[1]{\left|{#1}\right\rangle}$$
$$\newcommand{\bra}[1]{\left\langle{#1}\right|}$$
$$\newcommand{\braket}[2]{\left\langle{#1}\middle|{#2}\right\rangle}$$

# Quantum Decision Making

Decision-making models are central in economics and psychology. Finding suitable models which can approach human decision-making behavior is, as expected, a very challenging task. Economics has for a long time assumed models which are based on a particular axiomatic skeleton. The axioms are “reasonable” approximations of general decision situations.
A central starting point in preference modeling in economics consists in proving that there exists an equivalence between the preference relation of an object x over an object $y$, denoted as $x > y$, if and only if there exists a utility function $u(.)$, which maps a set of objects into $R$, such that $u(x) > u(y)$. This is not an easy task. However, as any good microeconomic theory textbook will show for a finite set of objects $X$ (to which $x$ and $y$ belong), the equivalence is relatively easy to show [[D01]](https://doi.org/10.1017/CBO9781139003261).
In psychology, Decision-Making is regarded as the cognitive process resulting in selecting a belief or a method of action among several possible alternative options. It could be either rational or irrational. Decision-making process is a reasoning process based on assumptions of the decision-maker's values, preferences, and beliefs. Every decision-making process produces a final choice, which may or may not prompt action [[D02]](https://en.wikipedia.org/wiki/Decision-making).
One famous example of decision-making theory is the prison's dilemma. One famous example of decision-making theory is the prison's dilemma. A prisoner's dilemma is a situation where individual decision-makers always have an incentive to choose to create a less than optimal outcome for the individuals as a group.

## Prison's Dilemma
The prisoner's dilemma is a standard example of a game analyzed in game theory that shows why two completely rational individuals might not cooperate, even if it appears that it is in their best interests to do so. It was originally framed by Merrill Flood and Melvin Dresher while working at RAND in 1950. Albert W. Tucker formalized the game with prison sentence rewards and named it "prisoner's dilemma" [[D03]](https://en.wikipedia.org/wiki/Prisoner%27s_dilemma).
The game is presented as:
Two members of a criminal organization are arrested and imprisoned. Each prisoner is in solitary confinement without communicating with the other. The prosecutors lack sufficient evidence to convict the pair on the principal charge, but they have enough to convict both on a lesser charge. Simultaneously, the prosecutors offer each prisoner a bargain. Each prisoner can either betray the other by testifying that the other committed the crime or cooperate with the other by remaining silent. The possible outcomes are:
* If A and B betray the other, they serve two years in prison,
* If A betrays B, but B remains silent, A will be set free, and B will serve three years in prison,
* If A remains silent but B betrays A, A will serve three years in prison, and B will be set free,
* If A and B both remain silent, both of them will serve only one year in prison (on the lesser charge).

It is implied that the prisoners will have no opportunity to reward or punish their partner other than the prison sentences they get and that their decision by itself will not affect their reputation in the future. As betraying a partner offers a greater reward than cooperating with them, all purely rational self-interested prisoners will betray the other, meaning the only possible outcome for two purely rational prisoners is for them to betray each other, even though cooperation would yield a greater reward.

This section will show how a Prison's dilemma's game can be modeled based on a quantum walk algorithm.

## Prison's Dilemma as a quantum walk model
Iterated bipartite quantum games can be implemented in the discrete-time quantum walk on the line. This section studies a discrete-time quantum walk on a line with two particles defined as two agents. Classically, random walks with $K$ particles are equivalent to $K$ independent single-particle random walks. In the quantum case, though, a walk with $K$ particles may contain quantum correlation, thus offering a resource unavailable in the classical scenario, introducing exciting features. Also, in the case of identical particles, we have to consider the effects of quantum statistics, giving an additional feature to quantum walks that can also be exploited.
We model the game based on two rational agents chosen from a restricted set of two-qubit unitary operations.
In the following, we present a quantum version of the Prisoner’s Dilemma in which both players use mixed strategies as a specific example [[D04]](https://arxiv.org/abs/quant-ph/0607143).

## Discrete-time quantum walk
The Hilbert space of a quantum walk on a line comprises two parts, $H = H_x \otimes H_c$.
$H_c$ is spanned by two orthonormal states $\ket{0_i}, \ket{1_i}$ as “coin” subspace. The spatial subspace, $H_x$, is traversed by the orthonormal set of position eigenstates, $|\ket{x_i}$, with $x ∈ Z$ labeling discrete sites on a line and we use this symbol to distinguish the position eigenstates from coin space. Repeated application $H_c$ generates the evolution. The evolution is generated by repeated application of a composite unitary operator $U$ which implements a coin operation, followed by a conditional shift in the walker's position.
We model the Prison's dilemma as a game where two walkers play, and two coins can model their decision in each step.
The quantum walk with two walkers A,B takes place in a Hilbert space $H_{AB} = H_A \otimes H_B$, where $H_A$ and $H_B$. After $N$ steps, a pure state characterized by a density operator $\rho_0 = \ket{\psi (0)}\bra{\psi (0)}$ evolves to $\rho_N = U^N \rho_0 U^{\dagger N}$ with [[D05]](http://www.phys-info.org/uploads/3/8/1/3/3813936/physreva.74.042304.pdf)

\begin{equation}
    U = S(I \otimes U_c)
\end{equation}

where $U_c$ is a unitary operation in coin subspace, $I$ is the identity and $S$ is a shift operation in $H_{AB}$.

## Model Prison's dilemma by quantum walk
We can model the game based on the quantum walk algorithm. In this model, two players are considered two walkers, and their decision is made by tossing a coin for each player. We can see a graphic example of Prison's dilemma in the following figure

We model each condition based on the following codes.

As the figure shows, the player's decision is considered coin states where each decision confessor or being silent is considered getting Heads or Tails. When the prison decides to remain silent or confess, it will be coded as $\ket{0}$ or $\ket{1}$ respectively.
On the other hand, the years of conviction are coded by position states. Since we have four states for each prison, we need two qubits for each one. So, we have two qubits for two coins presenting the decision of each player and four qubits in position space. At first, we consider the decision coins are Hadamard coins where are presented by $U(\pi/2, 0, \pi)$, and each player's state initially is $\ket{00}$.

Now, we are going to simulate the system:


text:
# A mathemactical model for calculating profiles of cosmogenic nuclides concentration vs. depth

In this notebook we show the equations of the model that we'll use to calculate profiles of nuclides concentration vs. depth. A more detailled description of the theoretical background behind this model can be found in, e.g., [Braucher et al., 2003][1], [Siame et al., 2004][2] and [Rixhon et al., 2011][3].

The contribution of a specific particule type $p$ to the concentration of a particular nuclide can be written as:

$$P_p,\Lambda_p \rightarrow c_p(z,\epsilon,t,\rho) = \frac{P_p} {\frac{\rho \, \epsilon}{\Lambda_p} + \lambda} \cdot \exp\left[{-\frac{\rho \, z}{\Lambda_p}}\right] \left[1 - \exp\left\{-t \left(\frac{\rho \, \epsilon}{\Lambda_p} + \lambda \right)\right\}\right]$$

where $z$ (cm) is the depth below the surface, $\epsilon$ (cm yr$^{-1}$) is the erosion rate, $t$ (yr) is the exposure time, $\rho$ (g cm$^{-3}$) is the soil density, $\lambda$ (yr$^{-1}$) is the radioactive decay constant, $P_p$ (atom g$^{-1}$yr$^{-1}$) is the production rate of the nuclide at the surface due to the particule, and $\Lambda_p$ (g cm$^{-2}$) is the effective apparent attenuation length for that particule type.

The nuclide concentration is the sum of individual contributions plus inheritance (i.e., a function of the concentration of the nuclide at the initiation of the exposure scenario). For example, for $^{10}Be$, we have:

$$
\begin{align}
C^{\mathrm{10Be}}(z,\epsilon,t,\rho) = & P_n,\Lambda_n & \rightarrow & c_n(z,\epsilon,t,\rho) & + \\
& P_{\mu \mathrm{Slow}},\Lambda_{\mu \mathrm{Slow}} & \rightarrow & c_{\mu \mathrm{Slow}}(z,\epsilon,t,\rho) & + \\
& P_{\mu \mathrm{Fast}},\Lambda_{\mu \mathrm{Fast}} & \rightarrow & c_{\mu \mathrm{Fast}}(z,\epsilon,t,\rho) & + \\
& C^{\mathrm{10Be}}_{t0} \cdot e^{-\lambda t} & &
\end{align}
$$

where $C^{\mathrm{10Be}}_{t0}$ (atoms g$^{-1}$) is the $^{10}Be$ concentration at the initiation of the exposure scenario, and $c_n(z,\epsilon,t,\rho)$, $c_{\mu \mathrm{Slow}}(z,\epsilon,t,\rho)$ and $c_{\mu \mathrm{Fast}}(z,\epsilon,t,\rho)$ are the contributions of neutrons, slow muons and fast muons, respectively. Inheritance is assumed to be equal for all sediment particles taken along the depth profile, which is a simplistic assumption though reasonable if the depth profile formed during a single geomorphic event that involved particle mixing.

Here below is a Python/Numpy vectorized implementation of this model. The code is saved in the file `models.py` so that is can be re-used in other notebooks.

[1]: http://dx.doi.org/10.1016/S0012-821X(03)00205-X
[2]: http://dx.doi.org/10.1016/S0012-821X(04)00061-5
[3]: http://dx.doi.org/10.1016/j.quageo.2010.11.001

```python
%%writefile models.py

"""
Models for calculating profiles of cosmogenic nuclide
concentration vs. depth

Author: B. Bovy
"""

import math

import numpy as np

def _particle_contrib(depth, erosion, exposure, density,
                      nuclide, particle):
    """
    Contribution of a specific particle type
    to the nuclide concentration.
    """
    return (
        particle['prod_rate'] / ((density * erosion /
                                  particle['damping_depth'])
                                 + nuclide['rdecay']) *
        np.exp(-1. * density * depth / particle['damping_depth']) *
        (1. - np.exp(-1. * exposure * ((erosion * density /
                                       particle['damping_depth'])
                                       + nuclide['rdecay'])))
    )

def C_nuclide(depth, erosion, exposure, density,
              inheritance, nuclide, particles):
    """
    Calculate the concentration(s) of a nuclide
    at given depth(s) (generic function).

    Parameters
    ----------
    depth : float or array_like
        depth(s) below the surface [cm]
    erosion : float or array_like
        erosion rate(s) [cm yr-1]
    exposure : float or array_like
        exposure time [yr]
    density : float or array_like
        soil density [g cm-3]
    inheritance : float or array_like
        concentration of the nuclide at the
        initiation of the exposure scenario
        [atoms g-1]
    nuclide : dict
        nuclide parameters
    particles : [dict, dict, ...]
        parameters related to each particule type
        that contribute to the nuclide production

    Returns
    -------
    float or array-like (broadcasted)
        the nuclide concentration(s) [atoms g-1]

    Notes
    -----
    if arrays are given for several arguments, they must
    have the same shape or must be at least be broadcastable.

    `nucleide` must have the following key(s):
        'rdecay': radioactive decay [yr-1]

    each item in `particles` must have the following keys:
        'prod_rate': surface production rate [atoms g-1 yr-1]
        'damping_depth': effective apparent attenuation depth
                         [g cm-2]

    """
    return (
        np.sum((_particle_contrib(depth, erosion,
                                  exposure, density,
                                  nuclide, p)
                for p in particles),
               axis=0) +
        inheritance * np.exp(-1. * nuclide['rdecay'] * exposure)
    )

def C_10Be(depth, erosion, exposure, density, inheritance,
           P_0=5.):
    """
    A model for profiles of 10Be concentration vs. depth.

    Notes
    -----
    The following parameters are used
    (Braucher et al. 2003)

    10Be radioactive decay: log(2) / 1.36e6

    Contribution of
    - neutrons
        - production rate: 0.9785 * P_0
        - damping depth: 160
    - slow muons
        - production rate: 0.015 * P_0
        - damping depth: 1500
    - fast muons
        - production rate: 0.0065 * P_0
        - damping depth: 5300

    See Also
    --------
    C_nuclide

    """
    # nuclide parameters
    berillium10 = {'rdecay': math.log(2.) / 1.36e6}

    # particles parameters
    neutrons = {'prod_rate': 0.9785 * P_0,
                'damping_depth': 160.}
    slow_muons = {'prod_rate': 0.015 * P_0,
                  'damping_depth': 1500.}
    fast_muons = {'prod_rate': 0.0065 * P_0,
                  'damping_depth': 5300.}

    return C_nuclide(depth, erosion, exposure,
                     density, inheritance,
                     berillium10,
                     [neutrons, slow_muons, fast_muons])

def C_26Al(depth, erosion, exposure, density, inheritance,
           P_0=35.):
    """
    A model for profiles of 26Al concentration vs. depth.

    Notes
    -----
    The following parameters are used
    (Braucher et al. 2003)

    26Al radioactive decay: log(2) / 0.72e6

    Contribution of
    - neutrons
        - production rate: 0.9785 * P_0
        - damping depth: 160
    - slow muons
        - production rate: 0.015 * P_0
        - damping depth: 1500
    - fast muons
        - production rate: 0.0065 * P_0
        - damping depth: 5300

    See Also
    --------
    C_nuclide

    """
    # nuclide parameters
    aluminium26 = {'rdecay': math.log(2.) / 0.72e6}

    # particles parameters
    neutrons = {'prod_rate': 0.9785 * P_0,
                'damping_depth': 160.}
    slow_muons = {'prod_rate': 0.015 * P_0,
                  'damping_depth': 1500.}
    fast_muons = {'prod_rate': 0.0065 * P_0,
                  'damping_depth': 5300.}

    return C_nuclide(depth, erosion, exposure,
                     density, inheritance,
                     aluminium26,
                     [neutrons, slow_muons, fast_muons])

```

    Overwriting models.py

## Model behavior

Interactive plotting of the concentration profile interactively to show the model sensitivity (IPython widget).

```python
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from IPython.html import widgets

import models

sns.set_context('notebook')
```

```python
depths = np.linspace(20, 500, 100)

def plot_profile_model(**kwargs):

    C10Be = models.C_10Be(depths,
                          kwargs['erosion_rate'],
                          kwargs['exposure_time'],
                          kwargs['soil_density'],
                          kwargs['inheritance'])

    fig, ax = plt.subplots()

    ax.plot(C10Be, -depths)
    plt.setp(ax,
             xlim=[0, 4e5], ylim=[-500, 0],
             xlabel='10Be concentration [atoms g-1]',
             ylabel='-1 * depth [cm]')

    return fig

```

```python
we = widgets.FloatSliderWidget(min=0., max=3e-3, step=5e-4, value=1e-3,
                               description="erosion rate [cm yr-1]")
wt = widgets.FloatSliderWidget(min=0., max=3e5, step=5e4, value=1e5,
                               description="exposure time [yr]")
wd = widgets.FloatSliderWidget(min=1.8, max=2.3, step=0.1, value=2.,
                               description="soil density [g cm-3]")
wi = widgets.FloatSliderWidget(min=0., max=1e5, step=2.5e4, value=0.5e5,
                               description="inheritance [atoms g-1]")

w = widgets.interact(plot_profile_model,
                     erosion_rate=we,
                     exposure_time=wt,
                     soil_density=wd,
                     inheritance=wi)
```

```python

```


text:
# Characterization of Systems in the Spectral Domain

*This Jupyter notebook is part of a [collection of notebooks](../index.ipynb) in the bachelors module Signals and Systems, Communications Engineering, Universität Rostock. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*

## The Transfer Function

Different time-domain representations are used to characterize linear time-invariant (LTI) systems. For instance, the constant coefficients of ordinary partial differential equations (ODEs), the [impulse response](../systems_time_domain/impulse_response.ipynb) and the [step response](../systems_time_domain/step_response.ipynb). The [transfer function](../systems_time_domain/eigenfunctions.ipynb#Transfer-Function) denotes the complex factor $H(s)$ an eigenfunction $e^{s t}$ at the input of an LTI system is weighted with when passing through the system. Since the Laplace transform decomposes signals with respect to these eigenfunctions, the transfer function constitutes a representation of an LTI system in the spectral (Laplace or Fourier) domain. The links between the transfer function and the time-domain representations of LTI systems are discussed in this section.

### Relation to Impulse Response and In-/Output Signal

It was already outlined in the context of the [convolution theorem of the Fourier](../fourier_transform/theorems.ipynb#Convolution-Theorem) and [Laplace transform](../laplace_transform/theorems.ipynb#Convolution-Theorem), that an LTI system can be represented equivalently in the temporal and spectral domain. In the temporal domain the output $y(t)$ of an LTI system is given by convolving the input signal $x(t)$ with the [impulse response](../systems_time_domain/impulse_response.ipynb) $h(t)$ of the system

\begin{equation}
y(t) = x(t) * h(t)
\end{equation}

Applying the convolution theorem of the Laplace transform yields

\begin{equation}
Y(s) = X(s) \cdot H(s)
\end{equation}

where $X(s) = \mathcal{L} \{ x(t) \}$ and $Y(s) = \mathcal{L} \{ y(t) \}$ are the Laplace transforms of the input and output signal, respectively. The transfer function $H(s)$ is given as the Laplace transform of the impulse response $h(t)$

\begin{equation}
H(s) = \mathcal{L} \{ h(t) \}
\end{equation}

This may also be concluded from the [previously derived link between the impulse response and the transfer function](../systems_time_domain/eigenfunctions.ipynb#Link-between-Transfer-Function-and-Impulse-Response). The derived results show that an LTI system can be fully characterized either in the temporal or spectral domain by its impulse response $h(t)$ or transfer function $H(s)$ or $H(j \omega)$, respectively.

It can furthermore be concluded that the transfer function is given as the quotient between output $Y(s)$ and input signal $X(s)$ in the Laplace domain

\begin{equation}
H(s) = \frac{Y(s)}{X(s)}
\end{equation}

for $X(s) \neq 0$. The same relations hold in the Fourier domain if the region of convergence (ROC) of the impulse response includes the imaginary axis $\Re \{ s \} = 0$. In general, the transfer function can be derived by divison of the spectra of the output and input signal. This can be used to measure the transfer function of a system by specific input signals that fulfill $X(s) \neq 0$ with $X(s) \approx$ constant.

**Example**

The [transfer function $H(s)$ of the 2nd order low-pass filter](../systems_time_domain/network_analysis.ipynb#Transfer-Function) was derived as

\begin{equation}
H(s) = \frac{1}{C L s^2 + C R s + 1}
\end{equation}

According to above findings, the impulse response $h(t)$ of the low-pass is given by inverse Laplace transform of the transfer function $H(s)$. First the transfer function is defined in `SymPy`

```python
import sympy as sym
%matplotlib inline
sym.init_printing()

R, L, C = sym.symbols('R L C', positive=True)
s = sym.symbols('s', complex=True)

H = 1/(C*L*s**2 + C*R*s + 1)
H
```

Then the inverse Laplace transform is computed for the specific normalized values $L = .5$, $R = 1$, $C = .4$ of the network elements

```python
t = sym.symbols('t', real=True)
RLC = {R: 1, L: sym.Rational('.5'), C: sym.Rational('.4')}

h = sym.inverse_laplace_transform(H.subs(RLC), s, t)
h
```

This result is plotted for illustration

```python
sym.plot(h, (t, 0, 10), xlabel='$t$', ylabel='h(t)')
```

    <sympy.plotting.plot.Plot at 0x10eafff98>

Note, the result is equal to the [impulse response derived by explicit solution of the underlying ODE](../systems_time_domain/network_analysis.ipynb#Impulse-Response).

### Relation to Differential Equations

Linear ODEs with constant coefficients describing the relation between input $x(t)$ and output $y(t)$ represent an LTI system. In this case a system is characterized by the constants $\alpha_n$ and $\beta_m$ of the ODE

\begin{equation}
\sum_{n=0}^{N} \alpha_n \frac{d^n}{d t^n} y(t) = \sum_{m=0}^{M} \beta_m \frac{d^m}{d t^m} x(t)
\end{equation}

where $\alpha_N \neq 0$ with $N$ denoting the degree of the ODE. Initial values are discarded in the following. Applying a two-sided Laplace transform to the left- and right-hand-side and using the [differentiation theorem](../laplace_transform/theorems.ipynb#Differentiation-Theorem) yields

\begin{equation}
Y(s) \sum_{n=0}^{N} \alpha_n s^n = X(s) \sum_{m=0}^{M} \beta_m s^m
\end{equation}

Rearranging this result derives the transfer function $H(s)$ of the system as

\begin{equation}
H(s) = \frac{Y(s)}{X(s)} = \frac{\sum_{m=0}^{M} \beta_m s^m}{\sum_{n=0}^{N} \alpha_n s^n}
\end{equation}

The transfer function of an LTI system whose in- and output relation can be characterized by an ODE is given as rational function in $s$. The coefficients of its numerator and denominator polynomial are given by the coefficients $\beta_m$ and $\alpha_n$ of the ODE. The conversion between both representations of an LTI system is therefore straightforward.

**Example**

The coefficients for the 2nd order low-pass filter can be determined [from its ODE](../systems_time_domain/network_analysis.ipynb#Differential-Equation) or from above transfer function as

\begin{align}
\beta_0 &= 1 & \alpha_0 &= 1 \\
& & \alpha_1 &= C R \\
& & \alpha_2 &= C L
\end{align}

From theses coefficients the transfer function $H(s)$ is constructed in `SymPy`

```python
R, L, C = sym.symbols('R L C', positive=True)
beta = 1
alpha = L*C, R*C, 1

H = sym.Poly(beta, s) / sym.Poly(alpha, s)
H
```

Many software packages for signal processing feature functionality for the computation of the transfer function from the coefficients $\beta_m$ and $\alpha_n$ of a system. For instance, the Python signal processing toolbox [`scipy.signal`](http://docs.scipy.org/doc/scipy/reference/signal.html) provides a function to numerically evaluate the transfer function $H(j \omega)$ at a number of angular frequencies $\omega_i$ for a given set of coefficients. This is illustrated in the following for the normalized values $R$, $L$ and $C$ used above.

```python
from scipy import signal

R = 1
L = 0.5
C = 0.4
beta = 1
alpha = L*C, R*C, 1

w, H = signal.freqresp((beta, alpha))
```

This results in the vector `w` of discrete angular frequencies $\omega_i$ and the vector `H` with the values of $H(j \omega)$ at these frequencies. The latter is plotted with the Python plotting library [`matplotlib`](http://matplotlib.org/).

```python
import matplotlib.pyplot as plt

plt.plot(w, abs(H))
plt.xlabel('$\omega$ in $s^{-1}$')
plt.ylabel('$|H(j \omega)|$')
plt.grid();
```

### Determination from Internal Structure

The transfer function $H(s)$ can be determined from the internal structure of an LTI system, if known. Either by derivation of the ODE characterizing the system, as shown above, or by direct analysis of the system in the Laplace/Fourier domain. For passive electrical networks this can be conducted by representing the signals and elements of the network by their complex impedances $Z(s)$ and applying Kirchhoff's circuit laws. This was illustrated for instance by the previous [Laplace domain analysis of a 2nd order low-pass filter](../laplace_transform/network_analysis.ipynb#Example:-Second-Order-Low-Pass-Filter). Similar procedures can be applied to other types of LTI systems, e.g. mechanical systems.

The transfer function is then given as the quotient between the output signal $Y(s)$ and the input signal $X(s)$ in the Laplace/Fourier domain.

**Copyright**

This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Sascha Spors, Continuous- and Discrete-Time Signals and Systems - Theory and Computational Examples*.


text:
## 3.4 Comparison between AWC and ATC DGDs

In this section we compare the pros and cons between AWC and ATC DGDs.

### 3.4.1 Convergence comparison

There are two major differences in convergence performance between AWC and ATC.

#### 3.4.1.1 Limiting Bias.

From Sec.3.2 and Sec.3.3, we know that

\begin{align}
\lim_{k\to\infty} \frac{1}{n}\sum_{i=1}^n \|x_i^{(k)} - x^\star\|^2 &= O\big( \frac{\alpha^2 \rho^2 b^2}{(1-\rho)^2} \big) \hspace{1cm} \mbox{(ATC-DGD)} \\
\lim_{k\to\infty} \frac{1}{n}\sum_{i=1}^n \|x_i^{(k)} - x^\star\|^2 &= O\big( \frac{\alpha^2 b^2}{(1-\rho)^2} \big) \hspace{1cm} \mbox{(AWC-DGD)}
\end{align}

It is observed that ATC-DGD has an imporved limiting bias given the same step-size $\alpha$ because $\rho < 1$.

#### 3.4.1.2 Step-size range

It is established in literatures \[Refs\] that ATC-DGD can converge if $\alpha \le \frac{2}{L}$ in which $L$ is the Lipschitz constant of each local gradient $\nabla f_i(x)$. However, it is established in \[Refs\] that AWC-DGD can converge if $\alpha \le \frac{1+\lambda_n(W)}{L}$ in which $\lambda_n(W) \in (-1,1)$ is the smallest eigenvalue of the combination matrix $W$. Apparently, AWC-DGD have a narrow step-size ranges to guarantee convergence.

### 3.4.2 Compare convergence performance with BlueFog

In this section we focus on the distributed least square problem and compare the performance of ATC and AWC with BlueFog. For the details of least square problem, see Sec.3.2 and 3.3.

#### 3.4.2.1 Set up BlueFog

In the following code, you should be able to see the id of your CPUs. We use 8 CPUs to conduct the following experiment.

```
import ipyparallel as ipp

rc = ipp.Client(profile="bluefog")
dview = rc[:]  # A DirectView of all engines
dview.block = True
rc.ids
```

    [0, 1, 2, 3, 4, 5, 6, 7]

Let each agent import necessary modules and then initialize BlueFog. You should be able to see the printed information like:

> \[stdout:0\] Hello, I am 1 among 8 processes
>
> ...

```
%%px
import numpy as np
import bluefog.torch as bf
import torch
from bluefog.common import topology_util
import networkx as nx
import matplotlib.pyplot as plt

%matplotlib inline

bf.init()
print(f"Hello, I am {bf.rank()} among {bf.size()} processes")
```

    [stdout:0] Hello, I am 7 among 8 processes
    [stdout:1] Hello, I am 6 among 8 processes
    [stdout:2] Hello, I am 4 among 8 processes
    [stdout:3] Hello, I am 2 among 8 processes
    [stdout:4] Hello, I am 1 among 8 processes
    [stdout:5] Hello, I am 5 among 8 processes
    [stdout:6] Hello, I am 3 among 8 processes
    [stdout:7] Hello, I am 0 among 8 processes

#### 3.4.2.2 Generate local data $A_i$ and $b_i$

```
%%px

def generate_data(m, n, x_o):

    A = torch.randn(m, n).to(torch.double)
    ns = 0.1 * torch.randn(m, 1).to(torch.double)
    b = A.mm(x_o) + ns

    return A, b
```

#### 3.4.2.3 Distributed gradient descent method

```
%%px

def distributed_grad_descent(A, b, maxite=5000, alpha=1e-1):

    x_opt = torch.zeros(n, 1, dtype=torch.double)

    for _ in range(maxite):
        # calculate local gradient
        grad_local = A.t().mm(A.mm(x_opt) - b)

        # global gradient
        grad = bf.allreduce(grad_local, name="gradient")

        # distributed gradient descent
        x_opt = x_opt - alpha * grad

    grad_local = A.t().mm(A.mm(x_opt) - b)
    grad = bf.allreduce(grad_local, name="gradient")  # global gradient

    # evaluate the convergence of distributed gradient descent
    # the norm of global gradient is expected to 0 (optimality condition)
    global_grad_norm = torch.norm(grad, p=2)
    if bf.rank() == 0:
        print(
            "[Distributed Grad Descent] Rank {}: global gradient norm: {}".format(
                bf.rank(), global_grad_norm
            )
        )

    return x_opt
```

In the following code we run distributed gradient descent to achieve the global solution $x^\star$ to the optimization problem. To validate whether $x^\star$ is optimal, it is enough to examine $\frac{1}{n}\sum_{i=1}^n \nabla f_i(x^\star) = 0$.

```
%%px

m, n = 20, 5
x_o = torch.randn(n, 1).to(torch.double)
x_o = bf.allreduce(x_o, name="x_o")  # synchronize the global reference point
A, b = generate_data(m, n, x_o)
x_opt = distributed_grad_descent(A, b, maxite=500, alpha=1e-2)
```

    [stdout:7] [Distributed Grad Descent] Rank 0: global gradient norm: 4.450207821573741e-15

#### 3.4.2.4 AWC-DGD

```
%%px

def AWC_DGD_one_step(x, x_opt, A, b, alpha=1e-2):

    # one-step ATC-DGD.
    # The combination weights have been determined by the associated combination matrix.

    grad_local = A.t().mm(A.mm(x) - b)  # compute local grad
    x_new = bf.neighbor_allreduce(x) - alpha * grad_local  # AWC update

    # the relative error: |x^k-x_gloval_average|/|x_gloval_average|
    rel_error = torch.norm(x_new - x_opt, p=2) / torch.norm(x_opt, p=2)

    return x_new, rel_error
```

#### 3.4.2.5 ATC-DGD

```
%%px

def ATC_DGD_one_step(x, x_opt, A, b, alpha=1e-2):

    # one-step ATC-DGD.
    # The combination weights have been determined by the associated combination matrix.

    grad_local = A.t().mm(A.mm(x) - b)  # compute local grad
    y = x - alpha * grad_local  # adapte
    x_new = bf.neighbor_allreduce(y)  # combination

    # the relative error: |x^k-x_gloval_average|/|x_gloval_average|
    rel_error = torch.norm(x_new - x_opt, p=2) / torch.norm(x_opt, p=2)

    return x_new, rel_error
```

#### 3.4.2.6 Limiting bias comparison between AWC and ATC

```
%%px

# Set topology as exponential-two topology.
G = topology_util.ExponentialTwoGraph(bf.size())
bf.set_topology(G)

maxite = 300

x_awc = torch.zeros(n, 1, dtype=torch.double)  # Initialize x
rel_error_awc = torch.zeros((maxite, 1))

x_atc = torch.zeros(n, 1, dtype=torch.double)  # Initialize x
rel_error_atc = torch.zeros((maxite, 1))

for ite in range(maxite):

    if bf.rank() == 0:
        if ite % 10 == 0:
            print("Progress {}/{}".format(ite, maxite))

    x_awc, rel_error_awc[ite] = AWC_DGD_one_step(
        x_awc, x_opt, A, b, alpha=2e-3
    )  # you can adjust alpha to different values
    x_atc, rel_error_atc[ite] = ATC_DGD_one_step(
        x_atc, x_opt, A, b, alpha=2e-3
    )  # you can adjust alpha to different values
```

    [stdout:7]
    Progress 0/300
    Progress 10/300
    Progress 20/300
    Progress 30/300
    Progress 40/300
    Progress 50/300
    Progress 60/300
    Progress 70/300
    Progress 80/300
    Progress 90/300
    Progress 100/300
    Progress 110/300
    Progress 120/300
    Progress 130/300
    Progress 140/300
    Progress 150/300
    Progress 160/300
    Progress 170/300
    Progress 180/300
    Progress 190/300
    Progress 200/300
    Progress 210/300
    Progress 220/300
    Progress 230/300
    Progress 240/300
    Progress 250/300
    Progress 260/300
    Progress 270/300
    Progress 280/300
    Progress 290/300

```
rel_error_awc_exp2 = dview.pull("rel_error_awc", block=True, targets=0)
rel_error_atc_exp2 = dview.pull("rel_error_atc", block=True, targets=0)
```

```
import matplotlib.pyplot as plt

%matplotlib inline

plt.semilogy(rel_error_awc_exp2)
plt.semilogy(rel_error_atc_exp2)

plt.legend(["AWC-DGD", "ATC-DGD"], fontsize=16)

plt.xlabel("Iteration", fontsize=16)
plt.ylabel("Relative error", fontsize=16)
```

With the same step-size 2e-3, it is observed that ATC-DGD has smaller limiting bias than AWC-DGD, which is consistent with the conclusion in Sec.3.4.1.1.

#### 3.4.2.7 Step-size stability range comparison

```
%%px

# Set topology as exponential-two topology.
G = topology_util.ExponentialTwoGraph(bf.size())
bf.set_topology(G)

maxite = 20

x_awc = torch.zeros(n, 1, dtype=torch.double)  # Initialize x
rel_error_awc = torch.zeros((maxite, 1))

x_atc = torch.zeros(n, 1, dtype=torch.double)  # Initialize x
rel_error_atc = torch.zeros((maxite, 1))

for ite in range(maxite):

    if bf.rank() == 0:
        if ite % 10 == 0:
            print("Progress {}/{}".format(ite, maxite))

    x_awc, rel_error_awc[ite] = AWC_DGD_one_step(
        x_awc, x_opt, A, b, alpha=3e-2
    )  # you can adjust alpha to different values
    x_atc, rel_error_atc[ite] = ATC_DGD_one_step(
        x_atc, x_opt, A, b, alpha=3e-2
    )  # you can adjust alpha to different values
```

    [stdout:7]
    Progress 0/20
    Progress 10/20

```
rel_error_awc_exp2 = dview.pull("rel_error_awc", block=True, targets=0)
rel_error_atc_exp2 = dview.pull("rel_error_atc", block=True, targets=0)
```

```
import matplotlib.pyplot as plt

%matplotlib inline

plt.semilogy(rel_error_awc_exp2)
plt.semilogy(rel_error_atc_exp2)

plt.legend(["AWC-DGD", "ATC-DGD"], fontsize=16)

plt.xlabel("Iteration", fontsize=16)
plt.ylabel("Relative error", fontsize=16)
```

With the same step-size 3e-2, it is observed that AWC diverges but ATC still converges, which is consistent with the conclusion in Sec.3.4.1.1.

### 3.4.2 Runtime comparison

While AWC has a little bit worse convergence performance than ATC, its AWC structure can overlap computation with its communication. Recall the AWC-DGD recursion:

$$\begin{align}
x^{(k+1)} = \underbrace{\sum_{i=1}^n w_{ij} x_j^{(k)}}_{\mathrm{comm.}} - \alpha \underbrace{\nabla f_i(x_i^{(k)})}_{\mathrm{compt.}} \quad \mbox{(AWC-DGD)}
\end{align}$$

It is observed that the computation of the gradient and the communication can be conducted in a parallel manner. In contrast, the ATC-DGD recursion

$$\begin{align}
x^{(k+1)} = \underbrace{\sum_{i=1}^n w_{ij} \Big(x_j^{(k)} - \alpha \underbrace{\nabla f_j(x_j^{(k)})}_{\mathrm{compt.}}\Big)}_{\mathrm{comm.}} \quad \mbox{(ATC-DGD)}
\end{align}$$

has to compute first and then communicate. Therefore, AWC-DGD can be implemented in a more efficient manner than ATC-DGD by overlapping communication with computation.

BlueFog supports the communication and computation overlapping for AWC-DGD. The utility ```bf.neighbor_allreduce_nonblocking(tensor)``` enables the non-blocking operation of neighbor allreduce. A typical way to overlap model communication and gradient computation is as follows:

```python
# one-step of the non-blocking AWC-DGD
x_handle = bf.neighbor_allreduce_nonblocking(x) # step 1: conduct neighbor-allreduce in a non-blocking way
grad_locl = ComputeGrad()                       # step 2: compute the gradient
x_new = bf.wait(x_handle) - alpha*grad_local    # step 3: finish the AWC-DGD update
```

Since neighbor-allreduce is conducted in a non-blcoking way, steps 1 and 2 are carried out in an independent and parallel manner. The above non-blocking code can save significant time than blocking AWC-DGD or ATC-DGD when either communication or computation is time-consuming.

The following example shows how fast the non-blocking AWC-DGD is compared to blocking AWC/ATC-DGD when solving the same problem.

```
%%px

def NonBlocking_AWC_DGD_one_step(x, x_opt, A, b, alpha=1e-2):

    # one-step NBK-AWC-DGD.

    x_handle = bf.neighbor_allreduce_nonblocking(x)
    grad_local = A.t().mm(A.mm(x) - b)  # compute local grad
    x_new = bf.wait(x_handle) - alpha * grad_local  # AWC update

    # the relative error: |x^k-x_gloval_average|/|x_gloval_average|
    rel_error = torch.norm(x_new - x_opt, p=2) / torch.norm(x_opt, p=2)

    return x_new, rel_error
```

```
%%px

m, n = 1000, 1000
x_o = torch.randn(n, 1).to(torch.double)
x_o = bf.allreduce(x_o, name="x_o")  # synchronize the global reference point
A, b = generate_data(m, n, x_o)
x_opt = distributed_grad_descent(A, b, maxite=300, alpha=1e-3)
```

    [stdout:7] [Distributed Grad Descent] Rank 0: global gradient norm: 1.6877124312783498e-12

```
%%px

import time

G = topology_util.ExponentialTwoGraph(
    bf.size()
)  # Set topology as exponential-two topology.
bf.set_topology(G)

maxite = 1000
alpha = 5e-5
rel_error_dict = {}

for method in ["ATC", "AWC", "NBK-AWC"]:
    if bf.rank() == 0:
        print("\nRunning {}:".format(method))

    x = torch.zeros(n, 1, dtype=torch.double).to(torch.double)  # Initialize x
    rel_error = torch.zeros((maxite, 1))
    time_list = []

    start = time.time()
    for ite in range(maxite):

        if bf.rank() == 0:
            if ite % 200 == 0:
                print("Progress {}/{}".format(ite, maxite))

        # you can adjust alpha to different values
        if method == "ATC":
            x, rel_error[ite] = ATC_DGD_one_step(x, x_opt, A, b, alpha=alpha)
        elif method == "NBK-AWC":
            x, rel_error[ite] = NonBlocking_AWC_DGD_one_step(
                x, x_opt, A, b, alpha=alpha
            )
        else:
            x, rel_error[ite] = AWC_DGD_one_step(x, x_opt, A, b, alpha=alpha)
        time_list.append(time.time() - start)

    rel_error_dict[method] = rel_error.numpy().reshape(-1)
    rel_error_dict[method + "_time"] = np.asarray(time_list)

    end = time.time()

    if bf.rank() == 0:
        print("{} finishes in {} seconds.".format(method, end - start))
```

    [stdout:7]
    Running ATC:
    Progress 0/1000
    Progress 200/1000
    Progress 400/1000
    Progress 600/1000
    Progress 800/1000
    ATC finishes in 23.32351303100586 seconds.

    Running AWC:
    Progress 0/1000
    Progress 200/1000
    Progress 400/1000
    Progress 600/1000
    Progress 800/1000
    AWC finishes in 23.93889021873474 seconds.

    Running NBK-AWC:
    Progress 0/1000
    Progress 200/1000
    Progress 400/1000
    Progress 600/1000
    Progress 800/1000
    NBK-AWC finishes in 21.54670763015747 seconds.

It is observed that Non-blocking AWC-DGD is the fastest by overlapping communication and computation.

max_stars_repo_name: Bluefog-Lib/bluefog-tutorial

max_issues_repo_name: Bluefog-Lib/bluefog-tutorial

max_forks_repo_name: Bluefog-Lib/bluefog-tutorial


text:
# Homework 2

**Instructions:** Complete the notebook below. Download the completed notebook in HTML format. Upload assignment using Canvas.

**Due:** Jan. 23 at **11am**.

## Exercise: NumPy Arrays

Follow the instructions in the following cells.

```python
# Import numpy
import numpy as np
```

```python
# Use the 'np.arange()' function to create a variable called 'numbers1' that stores the integers
# 1 through (and including) 10
numbers1 = np.arange(1,11)

# Print the value of 'numbers1'
print(numbers1)
```

    [ 1  2  3  4  5  6  7  8  9 10]

```python
# Use the 'np.arange()' function to create a variable called 'numbers2' that stores the numbers
# 0 through (and including) 1 with a step increment of 0.01
step =0.01
numbers2 = np.arange(0,1+step,step)

# Print the value of 'numbers2'
print(numbers2)
```

    [0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1  0.11 0.12 0.13
     0.14 0.15 0.16 0.17 0.18 0.19 0.2  0.21 0.22 0.23 0.24 0.25 0.26 0.27
     0.28 0.29 0.3  0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.4  0.41
     0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.5  0.51 0.52 0.53 0.54 0.55
     0.56 0.57 0.58 0.59 0.6  0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69
     0.7  0.71 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.8  0.81 0.82 0.83
     0.84 0.85 0.86 0.87 0.88 0.89 0.9  0.91 0.92 0.93 0.94 0.95 0.96 0.97
     0.98 0.99 1.  ]

```python
# Print the 5th value of 'numbers2'. (Remember that the index starts counting at 0)
print(numbers2[4])
```

    0.04

```python
# Print the last value of 'numbers2'.
print(numbers2[-1])
```

    1.0

```python
# Print the first 12 values of 'numbers2'.
print(numbers2[:12])
```

    [0.   0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1  0.11]

```python
# Print the last 12 values of 'numbers2'.
print(numbers2[-12:])
```

    [0.89 0.9  0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1.  ]

```python
# Use the 'np.zeros()' function to create a variable called 'zeros' that is an array of 20 zeros
zeros = np.zeros(20)

# Print the value of 'zeros'
print(zeros)
```

    [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

```python
# Change the second value of 'zeros' to 1 and print
zeros[1] = 1

# Print the value of 'zeros'
print(zeros)
```

    [0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

## Exercise: Random Numbers

Follow the instructions in the following cells.

```python
# Set the seed of NumPy's random number generator to 126
np.random.seed(126)

# Create a variable called 'epsilon' that is an array containing 25 draws from
# a normal distribution with mean 4 and standard deviation 2
epsilon = np.random.normal(loc=4,scale=2,size=25)

# Print the value of epsilon
print(epsilon)
```

    [ 2.67538919  6.80573472  2.82000233  4.20989426  4.16082178 -1.647651
      5.03812875  5.01159741  4.39492279  3.85680303  6.91478091  0.0854267
      7.02367775  3.94735082  2.48401748  7.02291824  6.25058482  5.57263063
      1.67321674  4.54716569  6.58060317  5.30662456  1.16478298  3.95415024
      3.56110405]

```python
# Print the mean of 'epsilon'
print(np.mean(epsilon))
```

    4.136587121846035

```python
# Print the standard deviation of 'epsilon'
print(np.std(epsilon))
```

    2.182193135563433

## Exercise: The Cobb-Douglas Production Function

The Cobb-Douglas production function can be written in per worker terms as :
  \begin{align}
  y & = A k^{\alpha},
  \end{align}
where $y$ denotes output per worker, $k$ denotes capital per worker, and $A$ denotes total factor productivity or technology.

### Part (a)

On a single axis: plot the Cobb-Douglas production for $A$ = 0.8, 1, 1.2, and 1.4 with $\alpha$ = 0.35 and $k$ ranging from 0 to 10. Each line should have a different color. Your plot must have a title and axis labels. The plot should also contain a legend that clearly indicates which line is associated with which value of $A$ and does not cover the plotted lines.

```python
# Import the pyplot module from Matplotlib as plt
import matplotlib.pyplot as plt

# Select the Matlplotlib style sheet to use (Optional)
plt.style.use('classic')

# Use the '%matplotlib inline' magic command to ensure that Matplotlib plots are displayed in the Notebook
%matplotlib inline
```

```python
# Set capital share (alpha)
alpha = 0.35

# Create an array of capital values
k = np.arange(0,10,0.001)

# Plot production function for each of the given values for A

def cobbDouglas(A,k,alpha):
    '''Returns the value of the production function A*k**alpha

    Args:
        A (float):          TFP or productivity value
        k (NumPy ndarray):  Array of capital values
        alpha (float):      Capital share in production

     Returns:
         NumPy series
     '''
    return A*k**alpha

for A in [0.8,1,1.2,1.4]:
    plt.plot(k,cobbDouglas(A,k,alpha),lw=3,alpha = 0.65,label='$A='+str(A)+'$')

# Add x- and y-axis labels
plt.xlabel('capital')
plt.ylabel('output')
# Add a title to the plot
plt.title('Cobb-Douglas production function')
# Create a legend
plt.legend(loc='lower right')
# Add a grid
plt.grid(True)
```

**Question**

1. *Briefly* explain in words how increasing $A$ affects the shape of the production function.

**Answer**

1. Increasing $A$ increases the height of the the production function at every value of capital except for 0 so increasing $A$ increases the steepness of the production function. <!-- answer -->

### Part (b)

On a single axis: plot the Cobb-Douglas production for $\alpha$ = 0.1, 0.2, 0.3, 0.4, and 0.5 with $A$ = 1 and $k$ ranging from 0 to 10. Each line should have a different color. Your plot must have a title and axis labels. The plot should also contain a legend that clearly indicates which line is associated with which value of $\alpha$ and does not cover the plotted lines.

```python
# Set TFP (A)
A = 1

# Plot production function for each of the given values for alpha
for alpha in [0.1,0.2,0.3,0.4,0.5]:
    plt.plot(k,cobbDouglas(A,k,alpha),lw=3,alpha = 0.65,label='$\\alpha='+str(alpha)+'$')

# Add x- and y-axis labels
plt.xlabel('capital')
plt.ylabel('output')
# Add a title to the plot
plt.title('Cobb-Douglas production function')
# Create a legend
plt.legend(ncol=3,loc='lower right')
# Add a grid
plt.grid(True)
```

**Question**

1. *Briefly* explain in words how increasing $\alpha$ affects the shape of the production function.

**Answer**

1. Increasing $\alpha$ reduces the curvature of the production function for capital between 0 and 1 and increases the steepness of the production function for capital greater than 1. <!-- answer -->

### Exercise: The Cardioid

The cardioid is a shape described by the parametric equations:

  \begin{align}
  x & = a(2\cos \theta - \cos 2\theta), \\
  y & = a(2\sin \theta - \sin 2\theta).
  \end{align}

Construct a well-labeled graph of the cardiod for $a=4$ and $\theta$ in $[0,2\pi]$. Your plot must have a title and axis labels.

```python
# Construct data for x and y
a = 4
theta = np.arange(0,2*np.pi,0.001)
x = a*(2*np.cos(theta) - np.cos(2*theta))
y = a*(2*np.sin(theta) - np.sin(2*theta))

# Plot y against x
plt.plot(x,y,lw=3,alpha=0.65)

# Create x-axis label
plt.xlabel('x')

# Create y-axis label
plt.ylabel('y')

# Create title for plot
plt.title('Cardioid')

# Add a grid to the plot
plt.grid()
```

## Exercise: Unconstrained optimization

Consider the quadratic function:

\begin{align}
f(x) & = -7x^2 + 930x + 30
\end{align}

You will use analytic (i.e., pencil and paper) and numerical methods to find the the value of $x$ that maximizes $f(x)$. Another name for $x$ that maximizes $f(x)$ is the *argument of the maximum* value $f(x)$.

### Part (a): Analytic solution

Use standard calculus methods to solve for the value of $x$ that maximizes $f(x)$ to **five decimal places**. Use your answer to complete the sentence in the next cell.

 The value of $x$ that maximizes $f(x)$ is:

### Part (b):  Numerical solution

In the cells below, you will use NumPy to try to compute the argument of the maximum of $f(x)$.

```python
# Use np.arange to create a variable called 'x' that is equal to the numbers 0 through 100
# with a spacing between numbers of 0.1
x = np.arange(0,100,0.1)

# Create a variable called 'f' that equals f(x) at each value of the array 'x' just defined
f = -7*x**2 + 930*x + 30

# Use np.argmax to create a variable called xstar equal to the value of 'x' that maximizes the function f(x).
xstar = x[np.argmax(f)]

# Print the value of xstar
print(xstar)
```

    66.4

```python
# Use np.arange to create a variable called 'x' that is equal to the numbers 0 through 100
# with a spacing between numbers of 0.001
x = np.arange(0,100,0.001)

# Create a variable called 'f' that equals f(x) at each value of the array 'x' just defined
f = -7*x**2 + 930*x + 30

# Use np.argmax to create a variable called xstar equal to the value of 'x' that maximizes the function f(x).
xstar = x[np.argmax(f)]

# Print the value of xstar
print(xstar)
```

    66.429

```python
# Use np.arange to create a variable called 'x' that is equal to the numbers 0 through *50*
# with a spacing between numbers of 0.001
x = np.arange(0,50,0.001)

# Create a variable called 'f' that equals f(x) at each value of the array 'x' just defined
f = -7*x**2 + 930*x + 30

# Use np.argmax to create a variable called xstar equal to the value of 'x' that maximizes the function f(x).
xstar = x[np.argmax(f)]

# Print the value of xstar
print(xstar)
```

    49.999

### Part (c): Evaluation

Provide answers to the follow questions in the next cell.

**Questions**

1. How did the choice of step size in the array `x` affect the accuracy of the computed results in the first two cells of Part (b)?

2. What do you think is the drawback to decreasing the stepsize in `x`?

3. In the previous cell, why did NumPy return value for `xstar` that is so different from the solution you derived in Part (a)?

**Answers**

1. Choosing a smaller step size creates a finer grid of points for which the function $f$ is calculated and therefore the accuracy is improved with a smaller step size. <!-- answer -->
2. A smaller step size means that the variable `x` will have more elements and so storing the variable and using it in computations will use more memory and will slow and/or crash the program. <!-- answer -->
3. The range of values in the variable `x` doesn't include $x^*$. <!-- answer -->

## Exercise: Utility Maximization

Recall the two good utility maximization problem from microeconomics. Let $x$ and $y$ denote the amount of two goods that a person consumes. The person receives utility from consumption given by:
  \begin{align}
  u(x,y) & = x^{\alpha}y^{\beta}
  \end{align}
The person has income $M$ to spend on the two goods and the price of the goods are $p_x$ and $p_y$. The consumer's budget constraint is:
  \begin{align}
  M & = p_x x + p_y y
  \end{align}
Suppose that $M = 100$, $\alpha=0.25$, $\beta=0.75$, $p_x = 1$. and $p_y = 0.5$. The consumer's problem is to maximize their utility subject to the budget constraint. While this problem can easily be solved by hand, we're going to use a computational approach. You can also solve the problem by hand to verify your solution.

### Part (a)

Use the budget constraint to solve for $y$ in terms of $x$, $p_x$, $p_y$, and $M$. Use the result to write the consumer's utility as a function of $x$ only. Create a variable called `x` equal to an array of values from 0 to 80 with step size equal to 0.001 and a variable called `utility` equal to the consumer's utility. Plot the consumer's utility against $x$.

```python
# Assign values to the constants alpha, beta, M, px, py
alpha = 0.25
beta = 0.75
M = 100
px=1
py=0.5

# Create an array of x values
x = np.arange(0,80,0.001)

# Create an array of utility values
utility = x**alpha*((M-px*x)/py)**beta

# Plot utility against x.
plt.plot(x,utility,lw=3,alpha=0.65)

# x- and y-axis labels
plt.xlabel('x')
plt.ylabel('utility')

# Title
plt.title('u[x,y(x)]')

# Add grid
plt.grid()
```

### Part (b)

The NumPy function `np.max()` returns the highest value in an array and `np.argmax()` returns the index of the highest value. Print the highest value and index of the highest value of `utility`.

```python
print('highest utility value         :',np.max(utility))
print('index of highest utility value:',np.argmax(utility))
```

    highest utility value         : 95.84146563694088
    index of highest utility value: 25000

### Part (c)

Use the index of the highest value of utility to find the value in `x` with the same index and store value in a new variable called `xstar`. Print the value of `xstar`.

```python
# Create variable 'xstar' equal to value in 'x' that maximizes utility
xstar = x[np.argmax(utility)]

# Print value of 'xstar'
print('xstar:',xstar)
```

    xstar: 25.0

### Part (d)

Use the budget constraint to the find the implied utility-maximizing vaue of $y$ and store this in a variable called `ystar`. Print `ystar`.

```python
# Create variable 'ystar' equal to value in 'y' that maximizes utility
ystar = M/py - px*xstar/py

# Print value of 'xstar'
print('ystar:',ystar)
```

    ystar: 150.0