text:
### Terminology and notation

#### Inputs
- $ s_t $ - state
- $ o_t $ - observation
- $ a_t $ - action

#### Outputs
- $ \pi_\theta(a_t|o_t) $ - policy
- $ \pi_\theta(a_t|s_t) $ - policy (fully observed)

### Mappings
- $ r(s, a) $: reward function

### Markov Decision Processes

#### Markov chain
- $ M = \{S, T\} $
- $ S - $ state space, states $ s \in S $ (discrete or continuous)
- $ T - $ transition operator, transition probability, dynamics function, $ p(s_{t+1}|s_t) $
    - let $ \mu_{t, i} = p(s_t = i) $ $ \vec{\mu_t} $ is a vector of probabilities
    - let $ T_{i, j} = p(s_{t+1}=i \lvert s_t = j) $, then $ \vec{\mu_{t+1}} = T\vec{\mu_t}$
    - T is a linear operator applied on the current vector of state probabilities to produce the next vector of state\                 probabilities. (2d matrix)

#### Markov property
- Given $ s_t $, $ s_{t+1} $ is independent of $ s_{t-1} $

- Markov chain by itself does not allow us to specify a decision making problem because there is no notion of actions.\
    (pioneered in 1950s.)

#### Markov Decision Process
- $ M = \{S, A, T, r\} $
- $ S - $ state space, states $ s \in S $ (discrete or continuous)
- $ A - $ action space, actions $ a \in A $ (discrete or continuous)
- $ T - $ transition operator, transition probability, dynamics function, $ p(s_{t+1}|s_t, a_t) $ (now a tensor, 3d matrix)
- $ r - $ reward function, $ r: S x A \rightarrow \mathbb{R} $

#### Parially observed Markov Decision Process
- $ M = \{S, \mathcal{A}, \mathcal{O}, T, \mathcal{E}, r\} $
- $ S - $ state space, states $ s \in S $ (discrete or continuous)
- $ \mathcal{A} - $ action space, actions $ a \in \mathcal{A} $ (discrete or continuous)
- $ \mathcal{O} - $ observation space, observations $ o \in \mathcal{O} $ (discrete or continuous)
- $ T - $ transition operator, transition probability, dynamics function, $ p(s_{t+1}|s_t, a_t) $ (now a tensor, 3d matrix)
- $ \mathcal{E} - $ emission probability $ p(o_t|s_t)$
- $ r - $ reward function, $ r: S x \mathcal{A} \rightarrow \mathbb{R} $

- We will be making decisions based on observations without access to the true states.

### The goal of reinforcement learning

#### To learn a policy
- $ \pi_{\theta}(a|s): s \rightarrow a $

#### Simulate environment dynamics to produce the next state
- $ p(s^{'}| s, a): s, a \rightarrow s^{'} $

#### Joint  probability distribution over states and actions (trajectory distribution) for a finite time horizon
- $ p_\theta(s_1, a_1, \dots, s_T, a_T )$ = $ p(s_1) \prod\limits_{t=1}^T \pi_\theta(a_t|s_t)p(s_{t+1}|s_t, a_t)$

#### Optimal policy
- $ \theta^{\star} = arg \max\limits_{\theta} E_{T \sim p_{\theta}(\tau) }[\sum\limits_t r(s_t, a_t)] $
- $ \tau $ stands for the trajectory. trajectories are sequences of states and actions.
- $ \theta $ represents the parameter of a neural network.
- expectation accounts for the stochasticity of the policy, transition probabilities, and the initial state distribution.
- expected value under the trajectory distribution of the sum of rewards.

#### Finite horizon case: state-action marginal using augmented  markov state
\begin{align}
\theta^{\star} = arg \max\limits_{\theta} E_{T \sim p_{\theta}(\tau) }\left[ \sum\limits_t r(s_t, a_t) \right]
\end{align}

\begin{eqnarray*}
\theta^{\star} = arg \max\limits_{\theta} \sum\limits_t E_{(s_t, a_t) \sim p_{\theta}(s_t, a_t) }\left[ r(s_t, a_t) \right]
\end{eqnarray*}

#### Infinite horizon case: stationary distribution
\begin{eqnarray*}
\theta^{\star} = arg \max\limits_{\theta} \sum\limits_{t=1}^T E_{(s_t, a_t) \sim p_{\theta}(s_t, a_t) }\left[ r(s_t, a_t) \right]
\end{eqnarray*}

- Assumptions
    - ergodicity
    - $\mu = T\mu, (T-I)\mu = 0$
    - $\mu$ is eigenvector of T with eigenvalue 1!

### Expectations and stochastic systems
- In reinforcement learning, we almost always care about increasing expected values (expectations).
- Why reinforcement learning algorithms can optimize non-smooth reward functions?
    - mainly because the expected values of non-smooth reward functions are smooth in $ \theta $. (How?)

### Reinforcement Learning Algorithms

#### The anatomy of a reinforcement learning algorithm
- Always contain three basic parts
    1. Generate samples (i.e. run the policy) - trial stage.\
        (sample trajectories from the trajectory distribution, typically the one induced by our policy)
    2. Fit a model/ estimate the return. (either model-based (dynamics) or implicit model (value-function))
    3. Improve the policy. (Change the policy to make it better.)
- Examples - Different incarnations of the second and third stage
    1. (Policy gradient)
        - Run policy $ \rightarrow $ Evaluate policy $ \rightarrow $ Improve policy
        - generate trajectories $ \rightarrow $ sum up the rewards $ \rightarrow $ try to make good trajectories be more likely.\
        calculate and apply gradients.
    2. RL by backprop (Model-based)

### Types of algorithms
1. Policy gradients: directly differentiate the above objective
2. Value-based: estimate value function or Q-function of the optimal policy (no explicit policy)
3. Actor-critic: estimate value function or Q-function of the current policy, use it to improve policy.
4. Model-based RL: estimate the transition model, and then
    - use it for planning (no explicit policy)
        - Trajectory optimization/optimal control (primarily in continuous spaces) - essentially backpropagation to optimize\
            actions
        - Discrete planning in discrete action spaces - e.g., Monte Carlo tree search
    - use it to improve a policy (backpropogate gradients)
        - Requires some tricks to make it work
    - use the model to learn a value function
        - Dynamic Programming
        - Generate simulated experience for model free learner.
    - something else

#### Rollout, sample the policy i.e. run the policy one step at a time.

### Tradeoffs between algorithms
(Why can't we just learn one RL algorithm that works for all tasks)\
- Different tradeoffs
    - Sample efficiency: How many samples do we need from orange box
    - Stability and ease of use: Multiple hyperparameters to select
- Different assumptions
    - Stochastic or deterministic?
    - Continuous or discrete?
    - Episodic or infinite horizon?
- Different things are easy or hard in different settings
    - Easier to represent the policy?
    - Easier to represent the model?

#### Sample Efficiency
- How many samples do we need to get a good policy?
- Most important question: Is the algorithm off policy?
    - Off policy: able to improve the policy without generating new samples from that policy
    - On policy: each time the policy is changed, even a little bit, we need to generate new samples

##### Why would we use a less efficient algorithm?
Wall clock time is not the same as efficiency!

#### Stability and ease of use
- Does it converge?
- And if it converges, to what?
- And does it converge every time?

Why is any of this even a question???
- Supervised learning: almost always gradient descent
- Reinforcement learning: often not gradient descent

#### Comparison: assumptions
- common assumption #1: full observability
    - generally assumed by value function fitting methods
    - can be mitigated by adding recurrence
- common assumption #2: episodic learning
    - Often assumed by pure policy gradient methods
    - assumed by some model-based RL methods
- common assumption #3: continuity or smoothness
    - Assumed by some continuous value function learning methods
    - Often assumed by some model-based RL methods

### Examples of algorithms
1. Value-based: estimate value function or Q-function of the optimal policy (no explicit policy)
    - Q-learning, DQN
    - Temporal difference learning
    - Fitted value iteration
2. Policy gradients: directly differentiate the above objective
    - REINFORCE
    - Natural policy gradient
    - Trust region policy optimization (TRPO)
    - Proximal Policy Optimization (PPO)
3. Actor-critic: estimate value function or Q-function of the current policy, use it to improve policy.
    - Asynchronous advantage actor-critic (A3C)
    - Soft actor-critic (SAC)
    - DDPg
4. Model-based RL: estimate the transition model, and the then
    - Dyna
    - Guided policy search

```python

```

max_stars_repo_name: mohankumarSriram/bob-the-robot

max_issues_repo_name: mohankumarSriram/bob-the-robot

max_forks_repo_name: mohankumarSriram/bob-the-robot


text:
# Example 1
- Using *Algo 1*, design estimator $\hat \pi(N)$ and compute $\hat \pi(10000)$

```
import numpy as np

def mcpi(N):
  n = 0
  for i in range(N):
    x = np.random.uniform(-1,1)
    y = np.random.uniform(-1,1)
    if (x**2 + y**2 < 1):
      n += 1
  return 4*n/N

mcpi(10000)
```

    3.1212

# Example 3
Given i.i.d $\{\alpha_i : i \in 1, 2, ..., N\}$, we use
$$\bar \alpha_N = \frac{1}{N} \sum_{i=1}^N \alpha_i$$
as its estimator of the mean $\mathbb E[\alpha_1]$ and use
$$\beta_N = \frac{1}{N} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2$$
as the estimator of $Var(\alpha_1)$. Suppose $\alpha_1 \in L^4$, then
- Prove $\beta_N$ is biased.
- Prove that $\beta_N$ is consistent in $L^2$.
- Can you propose an unbiased estimator?

__Proof:__
- To prove $\beta_N$ is biased, need to show that $\mathbb E[\beta_N] - Var(\alpha_1)$ does not equal to 0.

Because $\{\alpha_i\}_{i\in\mathbb{N}}$ is i.i.d., we can deduce that

\begin{equation}
\begin{aligned}
\mathbb E[\bar \alpha_N] &= \mathbb E[\frac{1}{N} \sum_{i=1}^N \alpha_i] \\
&= \frac{1}{N} \sum_{i=1}^N \mathbb E[\alpha_i] \\
&= \frac{1}{N} \cdot N \cdot \mathbb E[\alpha_1] \\
&= \mathbb E[\alpha_1], \\
Var(\bar \alpha_N) &= Var(\frac{1}{N} \sum_{i=1}^N \alpha_i) \\
&= \frac{1}{N^2} \sum_{i=1}^N Var(\alpha_i) \\
&= \frac{1}{N^2} \cdot N \cdot Var(\alpha_1) \\
&= \frac{1}{N} Var(\alpha_1)\\
\mathbb E[\bar \alpha_N^2] &= Var(\bar \alpha_N) + \mathbb E[\bar \alpha_N]^2 \\
&= \frac{1}{N} Var(\alpha_1) + \mathbb E[\alpha_1]^2 \\
\mathbb E[\alpha_1^2] &= \mathbb E[\alpha_1]^2 + Var(\alpha_1)
\end{aligned}
\end{equation}

Now observing $Bias(\beta_N)$, we have

\begin{equation}
\begin{aligned}
Bias(\beta_N) &= \mathbb E[\beta_N] - Var(\alpha_1) \\
&= \mathbb E[\frac{1}{N} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2] - Var(\alpha_1) \\
&= \frac{1}{N} \mathbb E[\sum_{i=1}^N (\alpha_i^2- 2 \alpha_i \bar \alpha_N + \bar \alpha_N^2)] - Var(\alpha_1) \\
&= \frac{1}{N} \mathbb E[\sum_{i=1}^N \alpha_i^2 - 2N \cdot \bar \alpha_N^2 + N \cdot \bar \alpha_N^2] - Var(\alpha_1) \\
&= \frac{1}{N} \sum_{i=1}^N \mathbb E[\alpha_i^2] - \frac{1}{N}\mathbb E[N \cdot \bar \alpha_N^2] - Var(\alpha_1) \\
&= \mathbb E[\alpha_1^2] - \mathbb E[\bar \alpha_N^2]- Var(\alpha_1) \\
&= - \frac{1}{N} Var(\alpha_1) \neq 0.\\
\end{aligned}
\end{equation}

Therefore, $\beta_N$ is biased.

- To prove that $\beta_N$ is consistent in $L^2$, we need to show that $MSE(\beta_N) \to 0$, when $N \to 0$, which means that
$$\mathbb E[(\beta_N - Var(\alpha_1))^2] \to 0.$$
Now we observe $\mathbb E[(\beta_N - Var(\alpha_1))^2]$, we have:

\begin{equation}
\begin{split}
E[(\beta_N - Var(\alpha_1))^2] &= \mathbb E[\beta_N^2 - 2\beta_N Var(\alpha_1) + Var(\alpha_1)^2] \\
&= \mathbb E[\beta_N^2] + \frac{2-N}{N} Var(\alpha_1)^2 \\
&= Var(\beta_N) + \mathbb E[\beta_N]^2 + \frac{2-N}{N} Var(\alpha_1)^2 \\
&= Var(\beta_N) + \frac{1}{N^2} Var(\alpha_1)^2 \\
\end{split}
\end{equation}

Now, we will use $T := \frac{(n-1)S^2}{\sigma^2} \sim \chi (n-1)$, where $S = \frac{1}{N-1} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2$ is the sample variance and $Var(\alpha_1) = \sigma^2$. Plus, we have $\mathbb E[T] = n-1, Var(T) = 2(n-1)$.

\begin{equation}
\begin{split}
Var(\beta_N) &= Var(\frac{N-1}{N} S) \\
&= \frac{(N-1)^2}{N^2} Var(S) \\
&= \frac{(N-1)^2}{N^2} \cdot \frac{2Var(\alpha_1)^2}{N-1} \\
&= \frac{2(N-1)}{N^2} Var(\alpha_1)^2 \\
E[(\beta_N - Var(\alpha_1))^2] &= \frac{2(N-1)}{N^2} Var(\alpha_1)^2 + \frac{1}{N^2} Var(\alpha_1)^2 \\
&= \frac{2N-1}{N^2} Var(\alpha_1)^2 \to 0, N \to \infty\\
\end{split}
\end{equation}

Therefore, $\beta_N$ is consitent in $L^2$.

- The unbiased estimator:
$$\hat \beta_N = \frac{1}{N-1} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2$$

\begin{equation}
\begin{aligned}
Bias(\hat \beta_N) &= \mathbb E[\hat \beta_N] - Var(\alpha_1) \\
&= \frac{1}{N-1} \sum_{i=1}^N \mathbb E[\alpha_i^2] - \frac{1}{N-1}\mathbb E[N \cdot \bar \alpha_N^2] - Var(\alpha_1) \\
&= \frac{N}{N-1} \bigg (\mathbb E[\alpha_1^2] - \mathbb E[\bar \alpha_N^2] \bigg )- Var(\alpha_1) \\
&= Var(\alpha_1) - Var(\alpha_1) = 0.\\
\end{aligned}
\end{equation}
which shows that $\hat \beta_N$ is unbiased estimator.

__Q.E.D.__

```
def beta(N,M):
  n = 0
  m = 0
  pi_list = []
  for i in range(N):
    pi_list.append(mcpi(M))
    n += pi_list[i]
  pi_bar = n/N
  for i in range(N):
    m += (pi_list[i] - pi_bar)**2
  beta = m/N
  return beta

beta(100,10000)
```

    0.00030025286400000054

```
import matplotlib.pyplot as plt

y = []
for i in range (5,11):
  y.append(beta(100,2**i))
x = [2**n for n in range(5,11)]
print(y)
plt.loglog(x,y)
```

```
def indi_func(a,sign1,b,sign2,x):       # sign1 and sign 2 mean whether we can take the endpoint respectively
  if (sign1 == 0):                      # 0 means that endpoint is not included
    if (sign2 == 0):
      if (a<x<b):
        f = 1
      else:
        f = 0
    elif (sign2 == 1):                  # 1 means that endpoint is included
      if (a<x<=b):
        f = 1
      else:
        f = 0
  elif (sign1 == 1):
    if (sign2 == 0):
      if (a<=x<b):
        f = 1
      else:
        f = 0
    elif (sign2 == 1):
      if (a<=x<=b):
        f = 1
      else:
        f = 0
  return f                              # Actually, there is no need to discuss whether endpoints are included or not
                                        # Since P(ept) = 0, x is uniform distribution
```

```
def mcintegral(N):
  s = 0
  for i in range(N):
    y = np.random.uniform(0,1)
    h = 100 * indi_func(0,0,0.01,1,y) + indi_func(0.01,0,1,0,y)
    s += h
  return s/N

mcintegral(10000000)
```

    1.9879606

```
def beta_prime(N,M):
  n = 0
  m = 0
  pi_list = []
  for i in range(N):
    pi_list.append(mcintegral(M))
    n += pi_list[i]
  pi_bar = n/N
  for i in range(N):
    m += (pi_list[i] - pi_bar)**2
  beta_prime = m/N
  return beta_prime

y = []
for i in range (5,11):
  y.append(beta_prime(100,2**i))
x = [2**n for n in range(5,11)]

plt.loglog(x,y)
```

ext: ipynb


text:
# Assignment 2

**Submission deadline: last lab session before or on Wednesday 25.10.17**

**Points: 13 + 1 bonus points**

## Downloading this notebook

This assignment is an Jupyter notebook. Download it by cloning https://github.com/janchorowski/nn_assignments. Follow the instructions in its README for instructions. Whenever possible, add your solutions to the notebook.

Please email us about any problems with it - we will try to correct them quickly. Also, please do not hesitate to use GitHub’s pull requests to send us corrections!

# Problem 1 [2p] Bayes' Theorem

Bayes' theorem allows to reason about conditional probabilities of causes and their effects:

\begin{equation}
p(A,B)=p(A|B)p(B)=p(B|A)p(A)
\end{equation}

\begin{equation}
p(A|B) = \frac{p(B|A)p(A)}{p(B)}
\end{equation}

Bayes' theorem allows us to reason about probabilities of causes, when
we observe their results.  Instead of directly answering the hard
question $p(\text{cause}|\text{result})$ we can instead separately
work out the marginal probabilities of causes $p(\text{cause})$ and
carefully study their effects $p(\text{effect}|\text{cause})$.

Solve the following using Bayes' theorem.

1. **[1p]** There are two boxes on the table: box \#1 holds two
  black balls and eight red ones, box \#2 holds 5 black ones and
  5 red ones. We pick a box at random (with equal probabilities),
  and then a ball from that box.
  1. What is the probability, that the
  ball came from box \#1 if we happened to pick a red ball?

1. **[1p]** The government has started a preventive program of
  mandatory tests for the Ebola virus. Mass testing method is
  imprecise, yielding 1% of false positives (healthy, but the test
  indicates the virus) and 1% of false negatives (
  having the virus but healthy according to test results).
  As Ebola is rather infrequent, lets assume that it occurs in
  one in a million people in Europe.
  1. What is the probability,
  that a random European, who has been tested positive for Ebola
  virus, is indeed a carrier?
  2. Suppose we have an additional information, that the person has just
  arrived from a country where one in a thousand people is a carrier.
  How much will be the increase in probability?
  3. How accurate should be the test, for a 80% probability of true
  positive in a European?

## 1A
$A$ - red ball had been picked up

$B$ - ball was picked up from box #1

$p(B|A)$ - probability, that the ball came from box #1 if we happened to pick a red ball

$p(B)$ - probability, that the ball came from box #1

$p(A)$ - probability, that we happened to pick a red ball

$p(A|B)$ - probability, that we've picked up ball from box #1

```python
pA=7/20.
pB=.5
pAB=.5*2/10
pBA=pAB*pB/pA
pBA
```

    0.14285714285714288

## 2

$A$ - Choosen European is a carrier

$B$ - Choosen European is not a carrier

$+$ - Test result is positive

$-$ - Test result is negative
\begin{equation}
p(A|+) = \frac{p(+|A)p(A)}{p(+|A)p(A)+p(+|B)p(B)}
\end{equation}

#### A
$p(-|A) = 0.01$

$p(+|A) = 1-p(-|A)= 0.99$

$p(+|B) = 0.01$

$p(A) = 1/10^6$

$p(B)=1-p(A)=0.999999$

```python
pPA=1-.01
pPB=1-pPA
pA=1/10.**6
pB=1-pA
pAP=pPA*pA/(pPA*pA+pPB*pB)
print pAP*100,"%"
```

    0.00989902989507 %

#### B
$p(A) = 1/10^3$

```python
pPA=1-.01
pPB=1-pPA
pA=1/10.**3 #now we know that one of thousand is a carrier
pB=1-pA
pAP=pPA*pA/(pPA*pA+pPB*pB)
print pAP*100,"%"
```

    9.01639344262 %

#### C
$p(A) = 1/10^6$

$p(+|A) = 80%$

```python
pPA=.8 # 80% true positive
pPB=1-pPA
pA=1/10.**6
pB=1-pA
pAP=pPA*pA/(pPA*pA+pPB*pB)
print pAP*100,"%"
```

    0.000399998800004 %

# Problem 2 [2p + 1b] Naive Bayes Classifier

The Bayes' theorem allows us to construct a classifier in which we
model how the data is generated. Here we will describe a
simple and popular example of such a classifier called the naive
Bayes classifier.  Despite its simplicity It is quite effective for
classification of text documents (e.g. as spam and non-spam).

Let a document be a sequence of words $D=W_1,W_2,\ldots,W_n$
We will model generation of text documents as a two-stage process.
First, document category $C_j$ is drawn at random with probability
$p(C_j)$, also called the *a priori* probability.
To define the class-conditional probability
$p(D|C_j)$, we will make a simplifying (naive)
assumption, that every word in the document is drawn independently at
random with probability $p(W_i|C)$:

\begin{equation*}
  p(D|C_j) = p(W_1,W_2,\ldots,W_n | C_j) \approx p(W_1|C_j)p(W_2|C_j)\ldots p(W_n|C_j).
\end{equation*}

To infer the class of a document we apply the Bayes theorem:
\begin{equation*}   p(C_j|D) = \frac{p(D|C_j)p(C_j)}{p(D)} = \frac{p(C_j)p(W_1|C_j)p(W_2|C_j)\ldots p(W_n|C_j)}{p(D)}.
\end{equation*}
Please note that since we assumed only a finite number of classes,
we can compute the term $p(D)$ by making sure that the *a
posteriori probabilities* $p(C_j|D)$ sum to $1$ over all classes.

In this exercise we will try to mimic the language-guessing feature
of [Google Translate](https://translate.google.com/), although
on a much smaller scale.  We are given an input which is a
lower-case sequence of characters (such as *"some people like
pineapple on their pizza"*), and we determine whether the
sequence's language is English, Polish or Spanish.
We will treat each character as a separate observation.
The numbers are taken from [Wikipedia article on letter frequency](https://en.wikipedia.org/wiki/Letter_frequency#Relative_frequencies_of_letters_in_other_languages). We display the first few rows:

```python
import pandas as pd
from StringIO import StringIO

wiki_table = u"""English|French|German|Spanish|Portuguese|Esperanto|Italian|Turkish|Swedish|Polish|Dutch|Danish|Icelandic|Finnish|Czech\na|8.167|7.636|6.516|11.525|14.634|12.117|11.745|12.920|9.383|10.503|7.486|6.025|10.110|12.217|8.421\nb|1.492|0.901|1.886|2.215|1.043|0.980|0.927|2.844|1.535|1.740|1.584|2.000|1.043|0.281|0.822\nc|2.782|3.260|2.732|4.019|3.882|0.776|4.501|1.463|1.486|3.895|1.242|0.565|0|0.281|0.740\nd|4.253|3.669|5.076|5.010|4.992|3.044|3.736|5.206|4.702|3.725|5.933|5.858|1.575|1.043|3.475\ne|12.702|14.715|16.396|12.181|12.570|8.995|11.792|9.912|10.149|7.352|18.91|15.453|6.418|7.968|7.562\nf|2.228|1.066|1.656|0.692|1.023|1.037|1.153|0.461|2.027|0.143|0.805|2.406|3.013|0.194|0.084\ng|2.015|0.866|3.009|1.768|1.303|1.171|1.644|1.253|2.862|1.731|3.403|4.077|4.241|0.392|0.092\nh|6.094|0.737|4.577|0.703|0.781|0.384|0.636|1.212|2.090|1.015|2.380|1.621|1.871|1.851|1.356\ni|6.966|7.529|6.550|6.247|6.186|10.012|10.143|9.600|5.817|8.328|6.499|6.000|7.578|10.817|6.073\nj|0.153|0.613|0.268|0.493|0.397|3.501|0.011|0.034|0.614|1.836|1.46|0.730|1.144|2.042|1.433\nk|0.772|0.049|1.417|0.011|0.015|4.163|0.009|5.683|3.140|2.753|2.248|3.395|3.314|4.973|2.894\nl|4.025|5.456|3.437|4.967|2.779|6.104|6.510|5.922|5.275|2.564|3.568|5.229|4.532|5.761|3.802\nm|2.406|2.968|2.534|3.157|4.738|2.994|2.512|3.752|3.471|2.515|2.213|3.237|4.041|3.202|2.446\nn|6.749|7.095|9.776|6.712|4.446|7.955|6.883|7.987|8.542|6.237|10.032|7.240|7.711|8.826|6.468\no|7.507|5.796|2.594|8.683|9.735|8.779|9.832|2.976|4.482|6.667|6.063|4.636|2.166|5.614|6.695\np|1.929|2.521|0.670|2.510|2.523|2.755|3.056|0.886|1.839|2.445|1.57|1.756|0.789|1.842|1.906\nq|0.095|1.362|0.018|0.877|1.204|0|0.505|0|0.020|0|0.009|0.007|0|0.013|0.001\nr|5.987|6.693|7.003|6.871|6.530|5.914|6.367|7.722|8.431|5.243|6.411|8.956|8.581|2.872|4.799\ns|6.327|7.948|7.270|7.977|6.805|6.092|4.981|3.014|6.590|5.224|3.73|5.805|5.630|7.862|5.212\nt|9.056|7.244|6.154|4.632|4.336|5.276|5.623|3.314|7.691|2.475|6.79|6.862|4.953|8.750|5.727\nu|2.758|6.311|4.166|2.927|3.639|3.183|3.011|3.235|1.919|2.062|1.99|1.979|4.562|5.008|2.160\nv|0.978|1.838|0.846|1.138|1.575|1.904|2.097|0.959|2.415|0.012|2.85|2.332|2.437|2.250|5.344\nw|2.360|0.074|1.921|0.017|0.037|0|0.033|0|0.142|5.813|1.52|0.069|0|0.094|0.016\nx|0.150|0.427|0.034|0.215|0.253|0|0.003|0|0.159|0.004|0.036|0.028|0.046|0.031|0.027\ny|1.974|0.128|0.039|1.008|0.006|0|0.020|3.336|0.708|3.206|0.035|0.698|0.900|1.745|1.043\nz|0.074|0.326|1.134|0.467|0.470|0.494|1.181|1.500|0.070|4.852|1.39|0.034|0|0.051|1.503\nà|0|0.486|0|0|0.072|0|0.635|0|0|0|0|0|0|0|0\nâ|0|0.051|0|0|0.562|0|0|0|0|0|0|0|0|0|0\ná|0|0|0|0.502|0.118|0|0|0|0|0|0|0|1.799|0|0.867\nå|0|0|0|0|0|0|0|0|1.338|0|0|1.190|0|0.003|0\nä|0|0|0.578|0|0|0|0|0|1.797|0|0|0|0|3.577|0\nã|0|0|0|0|0.733|0|0|0|0|0|0|0|0|0|0\ną|0|0|0|0|0|0|0|0|0|0.699|0|0|0|0|0\næ|0|0|0|0|0|0|0|0|0|0|0|0.872|0.867|0|0\nœ|0|0.018|0|0|0|0|0|0|0|0|0|0|0|0|0\nç|0|0.085|0|0|0.530|0|0|1.156|0|0|0|0|0|0|0\nĉ|0|0|0|0|0|0.657|0|0|0|0|0|0|0|0|0\nć|0|0|0|0|0|0|0|0|0|0.743|0|0|0|0|0\nč|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0.462\nď|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0.015\nð|0|0|0|0|0|0|0|0|0|0|0|0|4.393|0|0\nè|0|0.271|0|0|0|0|0.263|0|0|0|0|0|0|0|0\né|0|1.504|0|0.433|0.337|0|0|0|0|0|0|0|0.647|0|0.633\nê|0|0.218|0|0|0.450|0|0|0|0|0|0|0|0|0|0\në|0|0.008|0|0|0|0|0|0|0|0|0|0|0|0|0\nę|0|0|0|0|0|0|0|0|0|1.035|0|0|0|0|0\ně|0|0|0|0|0|0|0|0|0|0|0|0|0|0|1.222\nĝ|0|0|0|0|0|0.691|0|0|0|0|0|0|0|0|0\nğ|0|0|0|0|0|0|0|1.125|0|0|0|0|0|0|0\nĥ|0|0|0|0|0|0.022|0|0|0|0|0|0|0|0|0\nî|0|0.045|0|0|0|0|0|0|0|0|0|0|0|0|0\nì|0|0|0|0|0|0|0.030|0|0|0|0|0|0|0|0\ní|0|0|0|0.725|0.132|0|0|0|0|0|0|0|1.570|0|1.643\nï|0|0.005|0|0|0|0|0|0|0|0|0|0|0|0|0\nı|0|0|0|0|0|0|0|5.114|0|0|0|0|0|0|0\nĵ|0|0|0|0|0|0.055|0|0|0|0|0|0|0|0|0\nł|0|0|0|0|0|0|0|0|0|2.109|0|0|0|0|0\nñ|0|0|0|0.311|0|0|0|0|0|0|0|0|0|0|0\nń|0|0|0|0|0|0|0|0|0|0.362|0|0|0|0|0\nň|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0.007\nò|0|0|0|0|0|0|0.002|0|0|0|0|0|0|0|0\nö|0|0|0.443|0|0|0|0|0.777|1.305|0|0|0|0.777|0.444|0\nô|0|0.023|0|0|0.635|0|0|0|0|0|0|0|0|0|0\nó|0|0|0|0.827|0.296|0|0|0|0|1.141|0|0|0.994|0|0.024\nõ|0|0|0|0|0.040|0|0|0|0|0|0|0|0|0|0\nø|0|0|0|0|0|0|0|0|0|0|0|0.939|0|0|0\nř|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0.380\nŝ|0|0|0|0|0|0.385|0|0|0|0|0|0|0|0|0\nş|0|0|0|0|0|0|0|1.780|0|0|0|0|0|0|0\nś|0|0|0|0|0|0|0|0|0|0.814|0|0|0|0|0\nš|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0.688\nß|0|0|0.307|0|0|0|0|0|0|0|0|0|0|0|0\nť|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0.006\nþ|0|0|0|0|0|0|0|0|0|0|0|0|1.455|0|0\nù|0|0.058|0|0|0|0|0.166|0|0|0|0|0|0|0|0\nú|0|0|0|0.168|0.207|0|0|0|0|0|0|0|0.613|0|0.045\nû|0|0.060|0|0|0|0|0|0|0|0|0|0|0|0|0\nŭ|0|0|0|0|0|0.520|0|0|0|0|0|0|0|0|0\nü|0|0|0.995|0.012|0.026|0|0|1.854|0|0|0|0|0|0|0\nů|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0.204\ný|0|0|0|0|0|0|0|0|0|0|0|0|0.228|0|0.995\nź|0|0|0|0|0|0|0|0|0|0.078|0|0|0|0|0\nż|0|0|0|0|0|0|0|0|0|0.706|0|0|0|0|0\nž|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0.721"""
df = pd.read_table(StringIO(wiki_table), sep='|', index_col=0)
df.head()
```

      English
      French
      German
      Spanish
      Portuguese
      Esperanto
      Italian
      Turkish
      Swedish
      Polish
      Dutch
      Danish
      Icelandic
      Finnish
      Czech

      a
      8.167
      7.636
      6.516
      11.525
      14.634
      12.117
      11.745
      12.920
      9.383
      10.503
      7.486
      6.025
      10.110
      12.217
      8.421

      b
      1.492
      0.901
      1.886
      2.215
      1.043
      0.980
      0.927
      2.844
      1.535
      1.740
      1.584
      2.000
      1.043
      0.281
      0.822

      c
      2.782
      3.260
      2.732
      4.019
      3.882
      0.776
      4.501
      1.463
      1.486
      3.895
      1.242
      0.565
      0.000
      0.281
      0.740

      d
      4.253
      3.669
      5.076
      5.010
      4.992
      3.044
      3.736
      5.206
      4.702
      3.725
      5.933
      5.858
      1.575
      1.043
      3.475

      e
      12.702
      14.715
      16.396
      12.181
      12.570
      8.995
      11.792
      9.912
      10.149
      7.352
      18.910
      15.453
      6.418
      7.968
      7.562

Implement the language classifier and answer the following:

1. **[0.5p]** Naive Bayes can be implemented
    either by multiplying probabilities or by adding
    log-probabilities. Which one is better and why?
2. **[1.5p]** What is the language of the following phrases, according to the classifier (below in a code cell)? Assume equal prior language probabilities $P(C)$.
3. **[0-1 bonus]** What happens when a Naive Bayes classifier
      is applied to a document with out-of-vocabulary words? Propose
      some solutions. Relate them to the concept of Bayesian
      priors discussed during the lecture.

```python
sentences = [
    u"No dejes para mañana lo que puedas hacer hoy.",
    u"Przed wyruszeniem w drogę należy zebrać drużynę.",
    u"Żeby zrozumieć rekurencję, należy najpierw zrozumieć rekurencję.",
    u"Si vale la pena hacerlo vale la pena hacerlo bien.",
    u"Experience is what you get when you didn't get what you wanted.",
    u"Należy prowokować intelekt, nie intelektualistów.",
    u"Ĝojan Kristnaskon kaj feliĉan novan jaron.",
    u"Vos enfants sont très beaux. Ils sont adoptes?",
    u"Is het de rode of de groene knoop die de bom ontmantelen.",
    u"Se tem mesmo que cortar as unhas dos pés, faça o favor de deixar a cozinha.",
    u"Keine Elephanten in der Bar nach 8 Uhr!",
    u"Gostarias de ir a minha casa e fazer as coisas que de qualquer forma direi às pessoas que fizemos?",
    u"Vandaag heb ik 10 eetlepels soep gegeten",
    u"cuando tengo sed",
    u"hej boysy, looknijcie sobie przez windowsa"]
```

```python
#multiply version
langs = list(df)
letters = list(df.index)

# Normalize values to a probability distribution
df_norm = df.divide(df.values.sum(0) * len(langs))

from decimal import *
import operator
#set bigger precision
getcontext().prec = 28

def naive_bayes(sent, langs, df_norm):
    """Returns the most probable language of a sentence"""
    lang_score={}
    for lang in langs:
        score=Decimal(1)
        #chceck value of each letter
        for l in sent.lower():
            letter_score=1
            try:
                #we avoid all problems with non-letters characters
                letter_score = df_norm.loc[l.encode('utf-8'), lang]
            except:
                pass
            score*=Decimal(letter_score)
        lang_score[lang]=score
    #get last (biggest) value of sorted dict
    return sorted(lang_score.items(), key=operator.itemgetter(1))[-1][0]#+u" or "+ sorted(lang_score.items(), key=operator.itemgetter(1))[-2][0]

for sent in sentences:
    print naive_bayes(sent, langs, df_norm), ':', sent
```

    Spanish : No dejes para mañana lo que puedas hacer hoy.
    Polish : Przed wyruszeniem w drogę należy zebrać drużynę.
    Polish : Żeby zrozumieć rekurencję, należy najpierw zrozumieć rekurencję.
    Italian : Si vale la pena hacerlo vale la pena hacerlo bien.
    English : Experience is what you get when you didn't get what you wanted.
    Polish : Należy prowokować intelekt, nie intelektualistów.
    Esperanto : Ĝojan Kristnaskon kaj feliĉan novan jaron.
    French : Vos enfants sont très beaux. Ils sont adoptes?
    Dutch : Is het de rode of de groene knoop die de bom ontmantelen.
    Portuguese : Se tem mesmo que cortar as unhas dos pés, faça o favor de deixar a cozinha.
    German : Keine Elephanten in der Bar nach 8 Uhr!
    Portuguese : Gostarias de ir a minha casa e fazer as coisas que de qualquer forma direi às pessoas que fizemos?
    Dutch : Vandaag heb ik 10 eetlepels soep gegeten
    Spanish : cuando tengo sed
    Polish : hej boysy, looknijcie sobie przez windowsa

```python
#adding version
langs = list(df)
letters = list(df.index)

# Normalize values to a probability distribution
df_norm = df.divide(df.values.sum(0) * len(langs))

import math
import operator
getcontext().prec = 28

def naive_bayes(sent, langs, df_norm):
    """Returns the most probable language of a sentence"""
    lang_score={}
    for lang in langs:
        score=1.
        for l in sent.lower():
            letter_score=1.
            try:
                letter_score = df_norm.loc[l.encode('utf-8'), lang]
            except:
                pass
            if (letter_score != 0.):
                score+=math.log(letter_score)
            else:
                score = -float('Inf') #letter doesn't exist in language => sentence is not in this language
                break

        lang_score[lang]=score
    #print sorted(lang_score.items(), key=operator.itemgetter(1))
    return sorted(lang_score.items(), key=operator.itemgetter(1))[-1][0]

for sent in sentences:
    print naive_bayes(sent, langs, df_norm), ':', sent
```

    Spanish : No dejes para mañana lo que puedas hacer hoy.
    Polish : Przed wyruszeniem w drogę należy zebrać drużynę.
    Polish : Żeby zrozumieć rekurencję, należy najpierw zrozumieć rekurencję.
    Italian : Si vale la pena hacerlo vale la pena hacerlo bien.
    English : Experience is what you get when you didn't get what you wanted.
    Polish : Należy prowokować intelekt, nie intelektualistów.
    Esperanto : Ĝojan Kristnaskon kaj feliĉan novan jaron.
    French : Vos enfants sont très beaux. Ils sont adoptes?
    Dutch : Is het de rode of de groene knoop die de bom ontmantelen.
    Portuguese : Se tem mesmo que cortar as unhas dos pés, faça o favor de deixar a cozinha.
    German : Keine Elephanten in der Bar nach 8 Uhr!
    Portuguese : Gostarias de ir a minha casa e fazer as coisas que de qualquer forma direi às pessoas que fizemos?
    Dutch : Vandaag heb ik 10 eetlepels soep gegeten
    Spanish : cuando tengo sed
    Polish : hej boysy, looknijcie sobie przez windowsa

##### 3
It depend to frequency of this words. If they are not very common- it doesn't change classification. But if they are very common it's hard even to native speakers to recognized language of text.
To solve this we can use spell-checker which can cut off all non-vocabulary words. it shouldn't change 'value' of text much.

# Problem 3 [2p]
Given observations $x_1,\ldots,x_n$
  coming from a certain distribution,
  prove that MLE of a particular parameter of that distribution is equal to the sample mean $\frac{1}{n}\sum_{i=1}^n x_i$:
1. Bernoulli distribution with success probability $p$ and MLE $\hat{p}$,
2. Gaussian distribution $\mathcal{N}(\mu,\sigma)$ and MLE $\hat{\mu}$,
3. Poisson distribution $\mathit{Pois}(\lambda)$ and MLE $\hat{\lambda}$.

*NOTE: You can submit your solution on paper. To freeze the solution and spend late days, at minimum send a photo of your solution via e-mail.*

# Problem 4 [2p]

1. Find simple expressions for the following functions' derivatives with respect to vector $\mathbf{x}$.
    1. $\tanh(\mathbf{x})$
    2. $\sigma(\mathbf{x}) = \frac{1}{1 + e^{-\mathbf{x}}}$
2. Find the following functions' gradients with respect to vector $[x, y, z]^T$:
    1. $f_1([x, y, z]^T) = x + y$
    2. $f_2([x, y, z]^T) = xy$
    3. $f_3([x, y, z]^T) = x^2y^2$
    4. $f_4([x, y, z]^T) = (x + y)^2$
    5. $f_5([x, y, z]^T) = x^4 + x^2 y z + x y^2 z + z^4$
    6. $f_6([x, y, z]^T) = e^{x + 2y}$
    7. $f_7([x, y, z]^T) = \frac{1}{x y^2}$
    8. $f_8([x, y, z]^T) = ax + by + c$
    9. $f_9([x, y, z]^T) = \tanh(ax + by + c)$

# Problem 5 [2p]

Find the following functions' gradients or Jacobians with respect to vector $\mathbf{x}$, where $\mathbf{x}, \mathbf{b} \in \mathbb{R}^{n}$, $\mathbf{W} \in \mathbb{R}^{n \times n}$:

1. **[0.5p]** $\mathbf{W} \mathbf{x} + \mathbf{b}$
2. **[0.5p]** $\mathbf{x}^T \mathbf{W} \mathbf{x}$,
3. **[0.5p]** $\sigma(\mathbf{W} \mathbf{x} + \mathbf{b})$,
    with $\sigma$ applied element-wise. Hint: use the Hadamard product
    (denoted $\circ$) for element-wise matrix multiplication.
4. **[0.5p]** $-\log(S(\mathbf{x})_j)$, where $S$ is the
    softmax function
    (https://en.wikipedia.org/wiki/Softmax_function) and we are
    interested in the derivative over the $j$-th output of the
    Softmax.

# Problem 6 [2p] Linear Regression

1. Implement a function generating a dataset of $n$ points
  according to the following algorithm:
  1. Draw $n$ points $x \propto U(0;10)$ (uniformly distributed on $[0,10]$).
  2. Draw $n$ points $y \propto \mathcal{N}(1+20x-1.3x^2, 7)$
    (from a Gaussian distribution with $\mu=1+20x-1.3x^2$ and $\sigma=7$).

  Prepare a dataset of 30 elements and make a scatterplot of
  the expected value $y$ in function $x$.

2. Use linear regression to fit polynomials to the
  generated dataset. Fit polynomials of degrees zero (a constant line),
  one, two and three. An easy way to do it is to transform each data
  point $x$ into a vector of its powers $[1, x, x^2, \ldots, x^m]$.

  Plot the dataset and all regression curves on one figure.

  **Note:** The name _linear regression_ signifies that the
  hypothesis is linear with respect to parameters $\Theta$.
  However, the relationship between $x$ and $y$ is not constrained
  to a linear one. In this exercise it is a polynomial one.

```python
def make_dataset(N):
    X = TODO
    Y = TODO
    return X,Y

data = make_dataset(30)
scatter(TODO)
```

```python
from numpy.linalg import inv

def powers_of_X(X, degree):
    powers = np.arange(degree + 1).reshape(-1,1)
    return X**powers

def compute_polynomial(X, Theta):
    XP = powers_of_X(X, len(Theta) - 1) # len(Theta) x N
    Y = TODO
    return Y

plot_x_space = np.linspace(0,10,100)
scatter(data[0], data[1])
for degree in range(4):
    X = powers_of_X(data[0], degree) # Matrix d x N
    Y = data[1].reshape(1, -1)       # Matrix 1 x N
    Theta = TODO
    plot(plot_x_space, compute_polynomial(plot_x_space, Theta).ravel(),
         label="degree: %d" %(degree, ))
    print Theta, Theta.shape
legend(loc='lower right')
```

# Problem 7 [1p]

When the data set is small and highly dimensional (or when high degree polynomials are used) the linear regression solution may fit the noise in the data instead of capturing the general rule. We call this phenomenon overfitting and will discuss it in detail in a few lectures.

One way of preventing overfitting is to force the model's parameters to be small. We call this *regularization*. Consider the following cost function:

$$ J(\Theta) = \sum_{i=1}^N (y^{(i)} - \Theta^T x^{(i)})^2 + \frac{\lambda}{2} \Theta^T \Theta $$

Analyze datasets sampled using the following procedure:

1. $x \propto U(0;1)$: $x$ is sampled uniformly from the  $0-1$ range.
2. $y \propto \mathcal{N}(\mu=1+2x-5x^2 + 4x^3, \sigma=0.1)$: then
    $y$ is sampled from the Normal distribution with mean
    $\mu=1+2x-5x^2+4x^3$ and standard deviation $0.1$

Repeat 30 times an experiment in which you sample a new training
dataset, then fit polynomials of degree 0 to 14 and use $\lambda$
value from the set $\{0, 10^{-6}, 10^{-5}, 10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}\}$.

Plot the mean training and testing errors. What is the effect of increasing $\lambda$?

```python
def poly_fun(X, Theta):
    """Compute the value of polynomial with coefficients Theta for points in X"""
    #
    # TODO: Write body of the function.
    # Return a vector of values in points from X (the same shape as X).
    #

#
# The true polynomial relation:
# y(x) = 1 + 2x -5x^2 + 4x^3
#
# TODO: write down the proper coefficients
#
true_poly = np.array([1., 2., -5, 4])

def make_dataset(N, theta=true_poly, sigma=0.1):
    """ Sample a dataset """
    X = np.random.rand(N)
    Y = np.random.randn(N)*sigma + poly_fun(X, true_poly)
    return X,Y

train_data = make_dataset(30)
XX = np.linspace(0,1,100)
YY = poly_fun(XX, true_poly)
scatter(train_data[0], train_data[1], label='train data', color='r')
plot(XX, poly_fun(XX, true_poly), label='ground truth')
legend(loc='upper left')
```

```python
import sys
sys.path.append('..')

from common import plotting
```

```python
#please note: lambda is a reserved keyword in python, thus we use _lambda
def poly_fit(data, degree, _lambda):
    "Fit a polynomial of a given degree and weight decay parameter C"
    powers = np.arange(degree + 1.0).reshape(-1,1)
    X = data[0].reshape(1,-1)
    Y = data[1].reshape(1,-1)
    XX = X**powers
    #
    # TODO: implement the closed-form solution for Theta
    #
    # Please note that np.inv may be numerically unstable.
    # It is better to use np.linalg.solve or even a QR decomposition.
    #
    Theta = TODO
    return Theta
```

```python
lambdas = [0.0, 0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1]
degrees  = arange(15)

num_repetitions = 30
num_samples = 30

train_errors = np.zeros((len(lambdas), len(degrees)))
test_errors = np.zeros((len(lambdas), len(degrees)))

#sample a single dataset for all experiments
test_data = make_dataset(num_samples)

for repetition in xrange(num_repetitions):
    #sample a new training dataset for this repetition
    train_data = make_dataset(num_samples)
    #scatter(data[0], data[1])
    for degree_i, degree in enumerate(degrees):
        for lambda_i, _lambda in enumerate(lambdas):
            Theta = poly_fit(train_data, degree, _lambda)
            #
            # TODO: compute the mean training and test errors
            #
            train_errors[lambda_i, degree_i] += TODO
            test_errors[lambda_i, degree_i] += TODO
train_errors /=  num_repetitions
test_errors /= num_repetitions

figure(figsize=(10,20))
for lambda_i, _lambda in enumerate(lambdas):
    subplot(len(lambdas), 1, lambda_i+1)
    plot(degrees, train_errors[lambda_i,:], label='train')
    plot(degrees, test_errors[lambda_i,:],  label='test')
    ylim(0,0.1)
    title('lambda=%g'%(_lambda,))
```

max_stars_repo_name: iCarrrot/NN

max_issues_repo_name: iCarrrot/NN

max_forks_repo_name: iCarrrot/NN


text:
# 14 Linear Algebra: Singular Value Decomposition

One can always decompose a matrix $\mathsf{A}$

\begin{gather}
\mathsf{A} = \mathsf{U}\,\text{diag}(w_j)\,\mathsf{V}^{T}\\
\mathsf{U}^T \mathsf{U} = \mathsf{U} \mathsf{U}^T = 1\\
\mathsf{V}^T \mathsf{V} = \mathsf{V} \mathsf{V}^T = 1
\end{gather}

where $\mathsf{U}$ and $\mathsf{V}$ are orthogonal matrices and the $w_j$ are the _singular values_ that are assembled into a diagonal matrix $\mathsf{W}$.

$$
\mathsf{W} = \text{diag}(w_j)
$$

The inverse (if it exists) can be directly calculated from the SVD:

$$
\mathsf{A}^{-1} = \mathsf{V} \text{diag}(1/w_j) \mathsf{U}^T
$$

## Solving ill-conditioned coupled linear equations

```python
import numpy as np
```

### Non-singular matrix
Solve the linear system of equations

$$
\mathsf{A}\mathbf{x} = \mathbf{b}
$$

Using the standard linear solver in numpy:

```python
A = np.array([
        [1, 2, 3],
        [3, 2, 1],
        [-1, -2, -6],
    ])
b = np.array([0, 1, -1])
```

```python
np.linalg.solve(A, b)
```

    array([ 0.83333333, -0.91666667,  0.33333333])

Using the inverse from SVD:

$$
\mathbf{x} = \mathsf{A}^{-1} \mathbf{b}
$$

```python
U, w, VT = np.linalg.svd(A)
print(w)
```

    [7.74140616 2.96605874 0.52261473]

First check that the SVD really factors $\mathsf{A} = \mathsf{U}\,\text{diag}(w_j)\,\mathsf{V}^{T}$:

```python
U.dot(np.diag(w).dot(VT))
```

    array([[ 1.,  2.,  3.],
           [ 3.,  2.,  1.],
           [-1., -2., -6.]])

```python
np.allclose(A, U.dot(np.diag(w).dot(VT)))
```

    True

Now calculate the matrix inverse $\mathsf{A}^{-1} = \mathsf{V} \text{diag}(1/w_j) \mathsf{U}^T$:

```python
inv_w = 1/w
print(inv_w)
```

    [0.1291755  0.33714774 1.91345545]

```python
A_inv = VT.T.dot(np.diag(inv_w)).dot(U.T)
print(A_inv)
```

    [[-8.33333333e-01  5.00000000e-01 -3.33333333e-01]
     [ 1.41666667e+00 -2.50000000e-01  6.66666667e-01]
     [-3.33333333e-01 -1.08335035e-16 -3.33333333e-01]]

Check that this is the same that we get from `numpy.linalg.inv()`:

```python
np.allclose(A_inv, np.linalg.inv(A))
```

    True

Now, *finally* solve (and check against `numpy.linalg.solve()`):

```python
x = A_inv.dot(b)
print(x)
np.allclose(x, np.linalg.solve(A, b))
```

    [ 0.83333333 -0.91666667  0.33333333]

    True

```python
A.dot(x)
```

    array([-7.77156117e-16,  1.00000000e+00, -1.00000000e+00])

```python
np.allclose(A.dot(x), b)
```

    True

### Singular matrix

If the matrix $\mathsf{A}$ is *singular* (i.e., its rank (linearly independent rows or columns) is less than its dimension and hence the linear system of equation does not have a unique solution):

For example, the following matrix has the same row twice:

```python
C = np.array([
     [ 0.87119148,  0.9330127,  -0.9330127],
     [ 1.1160254,   0.04736717, -0.04736717],
     [ 1.1160254,   0.04736717, -0.04736717],
    ])
b1 = np.array([ 2.3674474,  -0.24813392, -0.24813392])
b2 = np.array([0, 1, 1])
```

```python
np.linalg.solve(C, b1)
```

NOTE: failure is not always that obvious: numerically, a matrix can be *almost* singular.

Try solving the linear system of equations

$$
\mathsf{D}\mathbf{x} = \mathbf{b}_1
$$
with matrix $\mathsf{D}$ below:

```python
D = C.copy()
D[2, :] = C[0] - 3*C[1]
D
```

    array([[ 0.87119148,  0.9330127 , -0.9330127 ],
           [ 1.1160254 ,  0.04736717, -0.04736717],
           [-2.47688472,  0.79091119, -0.79091119]])

```python
np.linalg.solve(D, b1)
```

    array([1.61493184e+00, 2.69013663e+16, 2.69013663e+16])

Note that some of the values are huge, and suspiciously like the inverse of machine precision? Sign of a nearly singular matrix.

**Note**: *Just because a function did not throw an exception it does not mean that the answer is correct.* **Always check your output!**

Now back to the example with $\mathsf{C}$:

#### SVD for singular matrices
If a matrix is *singular* or *near singular* then one can *still* apply SVD.

One can then compute the *pseudo inverse*

\begin{align}
\mathsf{A}^{-1} &= \mathsf{V} \text{diag}(\alpha_j) \mathsf{U}^T \\
\alpha_j &= \begin{cases}
 \frac{1}{w_j}, &\quad\text{if}\ w_j \neq 0\\
 0, &\quad\text{if}\ w_j = 0
 \end{cases}
\end{align}

i.e., any singular $w_j = 0$ is being "augmented" by setting

$$
\frac{1}{w_j} \rightarrow 0 \quad\text{if}\quad w_j = 0
$$

in $\text{diag}(1/w_j)$.

Perform the SVD for the singular matrix $\mathsf{C}$:

```python
U, w, VT = np.linalg.svd(C)
print(w)
```

    [1.99999999e+00 1.00000000e+00 2.46519033e-32]

Note the third value $w_2 \approx 0$: sign of a singular matrix.

Test that the SVD really decomposes $\mathsf{A} = \mathsf{U}\,\text{diag}(w_j)\,\mathsf{V}^{T}$:

```python
U.dot(np.diag(w).dot(VT))
```

    array([[ 0.87119148,  0.9330127 , -0.9330127 ],
           [ 1.1160254 ,  0.04736717, -0.04736717],
           [ 1.1160254 ,  0.04736717, -0.04736717]])

```python
np.allclose(C, U.dot(np.diag(w).dot(VT)))
```

    True

There are the **singular values** (let's say, $|w_i| < 10^{-12}$):

```python
singular_values = np.abs(w) < 1e-12
print(singular_values)
```

    [False False  True]

#### Pseudo-inverse

Calculate the **pseudo-inverse** from the SVD

\begin{align}
\mathsf{A}^{-1} &= \mathsf{V} \text{diag}(\alpha_j) \mathsf{U}^T \\
\alpha_j &= \begin{cases}
 \frac{1}{w_j}, &\quad\text{if}\ w_j \neq 0\\
 0, &\quad\text{if}\ w_j = 0
 \end{cases}
\end{align}

Augment:

```python
inv_w = 1/w
inv_w[singular_values] = 0
print(inv_w)
```

    [0.5 1.  0. ]

```python
C_inv = VT.T.dot(np.diag(inv_w)).dot(U.T)
print(C_inv)
```

    [[-0.04736717  0.46650635  0.46650635]
     [ 0.5580127  -0.21779787 -0.21779787]
     [-0.5580127   0.21779787  0.21779787]]

#### Solution for $\mathbf{b}_1$

Now solve the linear problem with SVD:

```python
x1 = C_inv.dot(b1)
print(x1)
```

    [-0.34365138  1.4291518  -1.4291518 ]

```python
C.dot(x1)
```

    array([ 2.3674474 , -0.24813392, -0.24813392])

```python
np.allclose(C.dot(x1), b1)
```

    True

Thus, using the pseudo-inverse $\mathsf{C}^{-1}$ we can obtain solutions to the equation

$$
\mathsf{C} \mathbf{x}_1 = \mathbf{b}_1
$$

However, $\mathbf{x}_1$ is not the only solution: there's a whole line of solutions that are formed by the special solution and a combination of the basis vectors in the *null space* of the matrix:

The (right) *kernel* or *null space* contains all vectors $\mathbf{x^0}$ for which

$$
\mathsf{C} \mathbf{x^0} = 0
$$

(The dimension of the null space corresponds to the number of singular values.) You can find a basis that spans the null space. Any linear combination of null space basis vectors will also end up in the null space when $\mathbf{A}$ is applied to it.

Specifically, if $\mathbf{x}_1$ is a special solution and $\lambda_1 \mathbf{x}^0_1 + \lambda_2 \mathbf{x}^0_2 + \dots$ is a vector in the null space then

$$
\mathbf{x} = \mathbf{x}_1 + ( \lambda_1 \mathbf{x}^0_1 + \lambda_2 \mathbf{x}^0_2 + \dots )
$$

is **also a solution** because

$$
\mathsf{C} \mathbf{x} = \mathsf{C} \mathbf{x_1} + \mathsf{C} ( \lambda_1 \mathbf{x}^0_1 + \lambda_2 \mathbf{x}^0_2 + \dots ) = \mathsf{C} \mathbf{x_1} + 0 = \mathbf{b}_1 + 0 = \mathbf{b}_1
$$

The $\lambda_i$ are arbitrary real numbers and hence there is an infinite number of solutions.

In SVD:

* The columns $U_{\cdot, i}$ of $\mathsf{U}$ (i.e. `U.T[i]` or `U[:, i]`) corresponding to non-zero $w_i$, i.e. $\{i : w_i \neq 0\}$, form the basis for the _range_ of the matrix $\mathsf{A}$.
* The columns $V_{\cdot, i}$ of $\mathsf{V}$ (i.e. `V.T[i]` or `V[:, i]`) corresponding to zero $w_i$, i.e. $\{i : w_i = 0\}$, form the basis for the _null space_ of the matrix $\mathsf{A}$.

```python
x1
```

    array([-0.34365138,  1.4291518 , -1.4291518 ])

The rank space comes from $\mathsf{U}^T$:

```python
U.T
```

    array([[-7.07106782e-01, -4.99999999e-01, -4.99999999e-01],
           [ 7.07106780e-01, -5.00000001e-01, -5.00000001e-01],
           [-2.47010760e-16, -7.07106781e-01,  7.07106781e-01]])

The basis vectors for the rank space (``~ bool_array`` applies a logical ``NOT`` operation to the entries in the boolean array so that we can pick out "not singular values"):

```python
U.T[~singular_values]
```

    array([[-0.70710678, -0.5       , -0.5       ],
           [ 0.70710678, -0.5       , -0.5       ]])

The null space comes from $\mathsf{V}^T$:

```python
VT
```

    array([[-0.8660254 , -0.35355339,  0.35355339],
           [-0.5       ,  0.61237244, -0.61237244],
           [-0.        , -0.70710678, -0.70710678]])

The basis vector for the null space:

```python
VT[singular_values]
```

    array([[-0.        , -0.70710678, -0.70710678]])

The component of $\mathbf{x}_1$ along the basis vector of the null space of $\mathsf{C}$ (here a 1D space) – note that this component is zero, i.e., the special solution lives in the rank space:

```python
x1.dot(VT[singular_values][0])
```

    2.220446049250313e-16

We can create a family of solutions by adding vectors in the null space to the special solution $\mathbf{x}_1$, e.g. $\lambda_1 = 2$:

```python
lambda_1 = 2
x1_1 = x1 + lambda_1 * VT[2]
print(x1_1)

np.allclose(C.dot(x1_1), b1)
```

    [-0.34365138  0.01493824 -2.84336536]

    True

Thus, **all** solutions are
```
x1 + lambda * VT[2]
```

#### Solution for  $\mathbf{b}_2$

The solution vector $x_2$ solves
$$
\mathsf{C}\mathbf{x}_2 = \mathbf{b}_2
$$

```python
b2
```

    array([0, 1, 1])

```python
x2 = C_inv.dot(b2)
print(x2)
print(C.dot(x2))
np.allclose(C.dot(x2), b2)
```

    [ 0.9330127  -0.43559574  0.43559574]
    [-4.4408921e-16  1.0000000e+00  1.0000000e+00]

    True

... and the general solution will again be obtained by adding any multiple of the null space basis vector.

#### Null space

The Null space is spanned by the following basis vectors (just one in this example):

```python
null_basis = VT[singular_values]
null_basis
```

    array([[-0.        , -0.70710678, -0.70710678]])

Show that
$$
\mathsf{C}\mathbf{x}^0 = 0
$$

```python
C.dot(null_basis.T)
```

    array([[ 0.0000000e+00],
           [-6.9388939e-18],
           [-6.9388939e-18]])

## SVD for fewer equations than unknowns
$N$ equations for $M$ unknowns with $N < M$:

* no unique solutions (underdetermined)
* $M-N$ dimensional family of solutions
* SVD: at least $M-N$ zero or negligible $w_j$: columns of $\mathsf{V}$ corresponding to singular $w_j$ span the solution space when added to a particular solution.

Same as the above [**Solving ill-conditioned coupled linear equations**](#Solving-ill-conditioned-coupled-linear-equations).

## SVD for more equations than unknowns
$N$ equations for $M$ unknowns with $N > M$:

* no exact solutions in general (overdetermined)
* but: SVD can provide best solution in the least-square sense
  $$
  \mathbf{x} = \mathsf{V}\, \text{diag}(1/w_j)\, \mathsf{U}^{T}\, \mathbf{b}
  $$
  where

  * $\mathbf{x}$ is a $M$-dimensional vector of the unknowns (parameters of the fit),
  * $\mathsf{V}$ is a $M \times N$ matrix
  * the $w_j$ form a square $M \times M$ matrix,
  * $\mathsf{U}$ is a $M \times N$ matrix (and $\mathsf{U}^T$ is a $N \times M$ matrix), and
  * $\mathbf{b}$ is the $N$-dimensional vector of the given values (data)

It can be shown that $\mathbf{x}$ minimizes the residual

$$
\mathbf{r} := |\mathsf{A}\mathbf{x} - \mathbf{b}|.
$$

where the matrix $\mathsf{A}$ will be described below and will contain the evaluation of the fit function for each data point in $\mathbf{b}$.

(For a $N \le M$, one can find $\mathbf{x}$ so that $\mathbf{r} = 0$ – see above.)

(In the following, we will switch notation and denote the vector of $M$ unknown parameters of the model as $\mathbf{a}$; this $\mathbf{a}$ corresponds to $\mathbf{x}$ above. $N$ is the number of observations.)

### Linear least-squares fitting

This is the *linear least-squares fitting problem*: Given $N$ data points $(x_i, y_i)$ (where $1 \le i \le N$), fit to a linear model $y(x)$, which can be any linear combination of $M$ functions of $x$.

For example, if we have $N$ functions $x^k$ with parameters $a_k$

$$
y(x) = a_1 + a_2 x + a_3 x^2 + \dots + a_M x^{M-1}
$$

or in general
$$
y(x) = \sum_{k=1}^M a_k X_k(x)
$$

The goal is to determine the $M$ coefficients $a_k$.

Define the **merit function**
$$
\chi^2 = \sum_{i=1}^N \left[ \frac{y_i - \sum_{k=1}^M a_k X_k(x_i)}{\sigma_i}\right]^2
$$
(sum of squared deviations, weighted with standard deviations $\sigma_i$ on the $y_i$).

Best parameters $a_k$ are the ones that *minimize $\chi^2$*.

*Design matrix* $\mathsf{A}$ ($N \times M$, $N \geq M$), vector of measurements $\mathbf{b}$ ($N$-dim) and parameter vector $\mathbf{a}$ ($M$-dim):

\begin{align}
A_{ij} &= \frac{X_j(x_i)}{\sigma_i}\\
b_i &= \frac{y_i}{\sigma_i}\\
\mathbf{a} &= (a_1, a_2, \dots, a_M)
\end{align}

The design matrix $\mathsf{A}$ contains the *predicted* values from the basis functions for all values $x_i$ of the independent variable $x$ for which we have measured data $y_i$.

Minimum occurs when the derivative vanishes:
$$
0 = \frac{\partial\chi^2}{\partial a_k} = \sum_{i=1}^N {\sigma_i}^{-2} \left[ y_i - \sum_{j=1}^M a_j X_j(x_i) \right] X_k(x_i), \quad 1 \leq k \leq M
$$
($M$ coupled equations)

To simplify the notation, define the $M \times M$ matrix
\begin{align}
\alpha_{kj} &= \sum_{i=1}^N \frac{X_k(x_i) X_j(x_i)}{\sigma_i^2}\\
\mathsf{\alpha} &= \mathsf{A}^T \mathsf{A}
\end{align}
and the vector of length $M$
\begin{align}
\beta_{k} &= \sum_{i=1}^N \frac{y_i X_k(x_i)}{\sigma_i^2}\\
\boldsymbol{\beta} &= \mathsf{A}^T \mathbf{b}
\end{align}

Then the $M$ coupled equations can be compactly written as
\begin{align}
\sum_{j=1}^{M} \alpha_{kj} a_j &= \beta_k\\
\mathsf{\alpha}\mathbf{a} = \boldsymbol{\beta}
\end{align}

$\mathsf{\alpha}$ and  $\boldsymbol{\beta}$ are known, so we have to solve this matrix equation for the vector of the unknown parameters $\mathbf{a}$.

#### Error estimates for the parameters

The inverse of $\mathsf{\alpha}$ is related to the uncertainties in the parameters:
$$
\mathsf{C} := \mathsf{\alpha}^{-1}
$$
in particular
$$
\sigma(a_i) = C_{ii}
$$
(and the $C_{ij}$ are the co-variances).

#### Solution of the linear least-squares fitting problem with SVD
We need to solve the overdetermined system of $M$ coupled equations
\begin{align}
\sum_{j=1}^{M} \alpha_{kj} a_j &= \beta_k\\
\mathsf{\alpha}\mathbf{a} = \boldsymbol{\beta}
\end{align}
We can solve the above equation with SVD.

SVD finds $\mathbf{a}$ that minimizes
$$
\chi^2 = |\mathsf{A}\mathbf{a} - \mathbf{b}|
$$
(proof in _Numerical Recipes_ Ch 2.) and so SVD is suitable to solve the equation above.

We can alternatively use SVD to directly solve the overdetermined system
$$
\mathsf{A}\mathbf{a} = \mathbf{b}
$$
and we should get the same answer.

The first approach might be preferrable because the matrix $\mathsf{\alpha}$ is a relatively small $M \times M$ matrix, i.e., its size depends on the number of parameters. The design matrix $\mathsf{A}$ is a $N \times M$ matrix and can be large (in one dimensions) because its size depends on the possibly large number $N$, the number of data points.

The errors are
$$
\sigma^2(a_j) = \sum_{i=1}^{M} \left(\frac{V_{ji}}{w_i}\right)^2
$$
(see also _Numerical Recipes_ Ch. 15) where $V_{ji}$ are elements of $\mathsf{V}$.

#### Example
Synthetic data

$$
y(x) = 3\sin x - 2\sin 3x + \sin 4x
$$

with noise $r$ added (uniform in range $-5 < r < 5$).

```python
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
matplotlib.style.use('ggplot')

import numpy as np
```

```python
def signal(x, noise=0):
    r = np.random.uniform(-noise, noise, len(x))
    return 3*np.sin(x) - 2*np.sin(3*x) + np.sin(4*x) + r
```

```python
X = np.linspace(-10, 10, 500)
Y = signal(X, noise=5)
```

```python
plt.plot(X, Y, 'r-', label="with noise")
plt.plot(X, signal(X, noise=0), 'k--', label="signal")
plt.xlabel("x")
plt.ylabel("y")
plt.legend(loc="best")
```

Define our fit function (the model) and the basis functions. We need the basis functions for setting up the problem and we will later use the fitfunction together with our parameter estimates to compare our fit to the true underlying function.

```python
def fitfunc(x, a):
    return a[0]*np.cos(x) + a[1]*np.sin(x) + \
           a[2]*np.cos(2*x) + a[3]*np.sin(2*x) + \
           a[4]*np.cos(3*x) + a[5]*np.sin(3*x) + \
           a[6]*np.cos(4*x) + a[7]*np.sin(4*x)

def basisfuncs(x):
    return np.array([np.cos(x), np.sin(x),
                     np.cos(2*x), np.sin(2*x),
                     np.cos(3*x), np.sin(3*x),
                     np.cos(4*x), np.sin(4*x)])
```

(Note that we could have used the `basisfuncs()` in `fitfunc()` – left as an exercise for the keen reader...)

Set up the $\mathsf{\alpha}$ matrix and the $\boldsymbol{\beta}$ vector (here we assume that all observations have the same error $\sigma = 1$):

```python
M = 8
sigma = 1.
alpha = np.zeros((M, M))
beta = np.zeros(M)
for x in X:
    Xk = basisfuncs(x)
    for k in range(M):
        for j in range(M):
            alpha[k, j] += Xk[k]*Xk[j]
for x, y in zip(X, Y):
    beta += y * basisfuncs(x)/sigma
```

Finally, solving the problem follows the same procedure as before:

Get the SVD:

```python
U, w, VT = np.linalg.svd(alpha)
V = VT.T
```

In this case, the singular values do not immediately show if any basis functions are superfluous (this would be the case for values close to 0).

```python
w
```

    array([296.92809624, 282.94804954, 243.7895787 , 235.7300808 ,
           235.15938555, 235.14838812, 235.14821093, 235.14821013])

... nevertheless, remember to routinely mask any singular values or close to singular values:

```python
w_inv = 1/w
w_inv[np.abs(w) < 1e-12] = 0
alpha_inv = V.dot(np.diag(w_inv)).dot(U.T)
```

Solve the system of equations with the pseudo-inverse:

```python
a_values = alpha_inv.dot(beta)
print(a_values)
```

    [ 0.22172712  3.19499849  0.18784697  0.00643614  0.08853002 -2.15177143
      0.07067379  0.87101436]

Compare the fitted values to the original parameters $a_j = 0, +3, 0, 0, 0, -2, 0, +1$.

The original parameters show up as 3.15, -2.04 and 1.08 but the other parameters also have appreciable values. Given that the noise was sizable, this is not unreasonable.

Compare the plot of the underlying true function ("signal", dashed line) to the model ("fit", solid line):

```python
plt.plot(X, fitfunc(X, a_values), 'b-', label="fit")
plt.plot(X, signal(X, noise=0), 'k--', label="signal")
plt.legend(loc="best", fontsize="small")
```

We get some spurious oscillations but overall the result looks reasonable.

```python

```

#### Direct calculation
Instead of solving the compact $M \times M$ matrix equation $\mathsf{\alpha}\mathbf{a} = \boldsymbol{\beta}$, we can try to directly solve the overdetermined $M \times N$ equation

$$
\mathsf{A}\mathbf{a} = \mathbf{b}
$$

Creating the design matrix is straight-forward (but because of the way that `basisfuncs()` returns values, we need to transpose the output to get the proper $N \times M$ matrix $\mathsf{A}$):

```python
A = np.transpose(basisfuncs(X))
b = Y
```

```python
A.shape
```

    (500, 8)

Calculate the pseudo-inverse $\mathsf{A}^{-1}$.

Note that we need to explicitly construct the matrix with the inverses of the singular values by filling a $M \times N$ matrix with $\text{diag}(w_i)$.

```python
U, w, VT = np.linalg.svd(A)
V = VT.T
singular_values = np.abs(w) < 1e-12
winv = 1/w
winv[singular_values] = 0
winvmat = np.zeros((V.shape[0], U.shape[0]))
winvmat[:len(winv), :len(winv)] = np.diag(winv)
Ainv = V.dot(winvmat).dot(U.T)
```

The singular values are all well behaved:

```python
w
```

    array([17.23160167, 16.8210597 , 15.61376248, 15.35350386, 15.33490742,
           15.33454884, 15.33454306, 15.33454304])

```python
V.shape, winvmat.shape, U.T.shape
```

    ((8, 8), (8, 500), (500, 500))

```python
Ainv.shape
```

    (8, 500)

```python
A.shape
```

    (500, 8)

Now solve directly
$$
\mathsf{A}^{-1} \mathbf{b} = \mathbf{a}
$$

```python
a = Ainv.dot(b)
a
```

    array([ 0.22172712,  3.19499849,  0.18784697,  0.00643614,  0.08853002,
           -2.15177143,  0.07067379,  0.87101436])

The parameter estimates are the same as above.

```python
a_values - a
```

    array([-9.15933995e-16,  2.22044605e-15, -6.38378239e-16, -3.44082401e-15,
           -1.36002321e-15, -1.33226763e-15,  4.16333634e-17,  3.55271368e-15])

and hence the plot looks the same:

```python
plt.plot(X, fitfunc(X, a_values), 'b-', label="fit")
plt.plot(X, signal(X, noise=0), 'k--', label="signal")
plt.legend(loc="best", fontsize="small")
```

```python

```


text:
```python
from sympy.physics.quantum.qubit import Qubit, IntQubit, IntQubitBra, matrix_to_qubit
from sympy.physics.quantum.dagger import Dagger
from sympy.physics.quantum.gate import Z
from sympy.physics.quantum.represent import represent
```

## Qubits

### We could have also called
Alternatively we could have also done
from sympy.physics.quantum.qubit import Qubit

```python
Qubit(0,0,0)
```

    |000>

```python
q = Qubit('0101')
q
```

    |0101>

```python
q.nqubits
```

    4

```python
len(q)
```

    4

```python
q.dimension
```

    4

```python
q.flip(1)
```

    |0111>

```python
q.qubit_values
```

    (0, 1, 0, 1)

```python
Dagger(q)
```

    <0101|

```python
type(Dagger(q))
```

    sympy.physics.quantum.qubit.QubitBra

```python
ip = Dagger(q)*q
ip
```

    <0101|0101>

```python
# Create a qubit for integer 5
q = IntQubit(Qubit('101'))
q
```

    |5>

```python
q.as_int()
```

    5

```python
q.nqubits
```

    3

```python
q.qubit_values
```

    (1, 0, 1)

```python
# Represent a state and then go back to its qubit form
q = Qubit('01')
matrix_to_qubit(represent(q))
```

    |01>

```python
# measure all, matrix to density
from sympy.physics.quantum.qubit import Qubit, measure_all, measure_partial, measure_all_oneshot, measure_partial_oneshot
from sympy.physics.quantum.gate import H, X, Y, Z
from sympy.physics.quantum.qapply import qapply
```

```python
c = H(0)*H(1)*Qubit('00')
c
```

    H(0)*H(1)*|00>

```python
q = qapply(c)
measure_all(q)
```

    [(|00>, 1/4), (|01>, 1/4), (|10>, 1/4), (|11>, 1/4)]

```python
# measure partial
measure_partial(q, (0,))
```

    [(sqrt(2)*|00>/2 + sqrt(2)*|10>/2, 1/2),
     (sqrt(2)*|01>/2 + sqrt(2)*|11>/2, 1/2)]

```python
# measure partial oneshot
measure_partial_oneshot(q, (0,))
```

    sqrt(2)*|01>/2 + sqrt(2)*|11>/2

```python
# measure all one shot
measure_all_oneshot(q)
```

    |10>

### Applying Hadamard gate to a single Qubit system

## Shor's Algorithm

From [Wikipedia](https://en.wikipedia.org/wiki/Shor%27s_algorithm), [**Also** read Peter Shor's arxiv research paper](https://arxiv.org/abs/quant-ph/9508027)

Shor's algorithm, named after mathematician Peter Shor, is a quantum algorithm (an algorithm that runs on a quantum computer) for integer factorization formulated in 1994. Informally it solves the following problem: given an integer N, find its prime factors.

On a quantum computer, to factor an integer N, Shor's algorithm runs in polynomial time (the time taken is polynomial in log N, which is the size of the input).[1] Specifically it takes quantum gates of order O((log N)2(log log N)(log log log N)) using fast multiplication,[2] demonstrating that the integer factorization problem can be efficiently solved on a quantum computer and is thus in the complexity class BQP. This is substantially faster than the most efficient known classical factoring algorithm, the general number field sieve, which works in sub-exponential time – about O(e1.9 (log N)1/3 (log log N)2/3).[3] The efficiency of Shor's algorithm is due to the efficiency of the quantum Fourier transform, and modular exponentiation by repeated squarings.

If a quantum computer with a sufficient number of qubits could operate without succumbing to noise and other quantum decoherence phenomena, Shor's algorithm could be used to break public-key cryptography schemes such as the widely used RSA scheme. RSA is based on the assumption that factoring large numbers is computationally intractable. So far as is known, this assumption is valid for classical (non-quantum) computers; no classical algorithm is known that can factor in polynomial time. However, Shor's algorithm shows that factoring is efficient on an ideal quantum computer, so it may be feasible to defeat RSA by constructing a large quantum computer. It was also a powerful motivator for the design and construction of quantum computers and for the study of new quantum computer algorithms. It has also facilitated research on new cryptosystems that are secure from quantum computers, collectively called post-quantum cryptography.

In 2001, Shor's algorithm was demonstrated by a group at IBM, who factored 15 into 3 × 5, using an NMR implementation of a quantum computer with 7 qubits.[4] After IBM's implementation, two independent groups implemented Shor's algorithm using photonic qubits, emphasizing that multi-qubit entanglement was observed when running the Shor's algorithm circuits.[5][6] In 2012, the factorization of 15 was performed with solid-state qubits.[7] Also in 2012, the factorization of 21 was achieved, setting the record for the largest number factored with Shor's algorithm.[8] In April 2012, the factorization of 143 was achieved, although this used adiabatic quantum computation rather than Shor's algorithm.[9] In November 2014, it was discovered that this 2012 adiabatic quantum computation had also factored larger numbers, the largest being 56153.[10][11]

### Plain english?

Let's say you have two prime numers. $p = 17$ and $q = 19$, then multiply these as follows as $pq = 323$.
Now lets subtract these two by 1 and multiply that as in:

$$ (p-1)*(q-1) = 16*18 = 288$$.

```python
p = 17
q = 19
A = p*q
B = (p-1)*(q-1)
print("Product of p and q is :",A)
print("Product of p-1 and q-1 is :", B)
```

    Product of p and q is : 323
    Product of p-1 and q-1 is : 288

#### What next?
Now, in a range of 1 and B(288), pick any number such that its LCD (lowest common denimonator) with A(323) is 1.
Let's say we pick 9?

What you'll realize is that you have to go a whole lot before you can find that number which is divisible by 9.

```python
num = 0
for x in range(1, 289):
    if (num % x) ==0:
        print ("Fact :")
```

```python
! python shor.py
```
	projectq
	--------
	Implementation of Shor's algorithm.
	Number to factor: 15

	Factoring N = 15: 00000001

	Bad luck: Found 15 and 1
Allocate : 109
CCR(0.19634954084936207) : 16
CCR(0.39269908169872414) : 32
CCR(0.7853981633974483) : 64
CCR(1.1780972450961724) : 8
CCR(1.3744467859455345) : 16
CCR(1.5707963267948966) : 104
CCR(1.9634954084936207) : 8
CCR(2.1598449493429825) : 16
CCR(2.356194490192345) : 24
CCR(2.552544031041707) : 16
CCR(2.748893571891069) : 32
CCR(3.141592653589793) : 192
CCR(3.5342917352885173) : 16
CCR(3.730641276137879) : 8
CCR(3.9269908169872414) : 24
CCR(4.123340357836604) : 8
CCR(4.319689898685965) : 16
CCR(4.71238898038469) : 88
CCR(4.908738521234052) : 8
CCR(5.105088062083414) : 16
CCR(5.497787143782138) : 56
CCR(5.890486225480862) : 16
CCR(6.086835766330224) : 8
CR(0.19634954084936207) : 144
CR(0.39269908169872414) : 288
CR(0.7853981633974483) : 432
CR(1.5707963267948966) : 568
CR(10.995574287564276) : 568
CR(11.780972450961723) : 432
CR(12.173671532660448) : 288
CR(12.370021073509811) : 144
CR(2.945243112740431) : 64
CR(3.141592653589793) : 64
CR(4.71238898038469) : 64
CR(5.497787143782138) : 64
CR(5.890486225480862) : 64
CSwap : 32
CX : 128
Deallocate : 109
H : 1408
Measure : 9
R(0.39269908169872414) : 64
R(0.7853981633974483) : 64
R(1.5707963267948966) : 64
R(3.141592653589793) : 64
R(3.3379421944391554) : 64
X : 130

Max. width (number of qubits) : 11.	projectq
	--------
	Implementation of Shor's algorithm.
	Number to factor: 21

	Factoring N = 21: 0000000000

	Factors found :-) : 3 * 7 = 21
Allocate : 166
CCR(0.09817477042468103) : 20
CCR(0.19634954084936207) : 40
CCR(0.39269908169872414) : 76
CCR(0.4908738521234052) : 20
CCR(0.7853981633974483) : 96
CCR(0.9817477042468103) : 22
CCR(1.0799224746714913) : 2
CCR(1.1780972450961724) : 30
CCR(1.2762720155208536) : 20
CCR(1.5707963267948966) : 189
CCR(1.6689710972195777) : 18
CCR(1.8653206380689396) : 20
CCR(1.9634954084936207) : 43
CCR(2.1598449493429825) : 2
CCR(2.356194490192345) : 53
CCR(2.552544031041707) : 30
CCR(2.945243112740431) : 9
CCR(3.141592653589793) : 300
CCR(3.3379421944391554) : 18
CCR(3.730641276137879) : 30
CCR(3.9269908169872414) : 73
CCR(4.123340357836604) : 1
CCR(4.319689898685965) : 23
CCR(4.417864669110646) : 10
CCR(4.614214209960008) : 9
CCR(4.71238898038469) : 111
CCR(5.006913291658733) : 10
CCR(5.105088062083414) : 30
CCR(5.203262832508095) : 1
CCR(5.3014376029327765) : 11
CCR(5.497787143782138) : 48
CCR(5.792311455056181) : 10
CCR(5.890486225480862) : 38
CCR(6.086835766330224) : 20
CCR(6.185010536754905) : 10
CR(0.09817477042468103) : 220
CR(0.19634954084936207) : 440
CR(0.39269908169872414) : 660
CR(0.7853981633974483) : 880
CR(1.5707963267948966) : 1191
CR(1.9634954084936207) : 100
CR(10.995574287564276) : 1091
CR(11.780972450961723) : 880
CR(12.173671532660448) : 660
CR(12.370021073509811) : 440
CR(12.468195843934492) : 220
CR(2.061670178918302) : 100
CR(3.141592653589793) : 100
CR(3.9269908169872414) : 100
CR(4.123340357836604) : 100
CSwap : 50
CX : 200
Deallocate : 166
H : 2602
Measure : 11
R(2.1598449493429825) : 100
R(2.356194490192345) : 100
R(3.141592653589793) : 100
R(4.221515128261284) : 100
R(4.319689898685965) : 100
R(4.71238898038469) : 100
X : 201

Max. width (number of qubits) : 13.

```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def Hadamard(n):
    def Hn(H=np.array([[1, 1], [1, -1]], dtype=np.complex64), n=n):
        if n > 1:
            return Hn(H=np.kron(np.array([[1, 1], [1, -1]], dtype=np.complex64), H), n=n-1)
        return H

    return Hn(n=n)

def QFT(t):
    Q = np.zeros(shape=(2 ** t, 2 ** t), dtype=np.complex64)
    N = 2 ** t
    for i in range(N):
        for j in range(N):
            Q[i][j] = np.exp(np.pi * 2j * ((i * j) % N) / N)

    return Q

N = 21
t = 9
H = Hadamard(t)

reg1 = np.zeros(shape=(2 ** t), dtype=np.complex64)
reg2 = np.ones(shape=(2 ** t), dtype=np.complex64)
reg1[0] = 1
reg1 = H.dot(reg1)

for i in range(2 ** t):
    reg2[i] = 2 ** i % N

r = reg2[0]

for i in range(2 ** t):
    if reg2[i] != r:
        reg1[i] = 0

Q = QFT(9)
reg1 = np.linalg.inv(Q).dot(reg1)
print (abs(reg1))
print (abs(reg1[0]))
print (abs(reg1[85]))
print (abs(reg1[86]))

fig, ax = plt.subplots( nrows=1, ncols=1 )
ax.plot(abs(reg1))
fig.savefig('plot.png')
plt.show()
```

```python

```

max_stars_repo_path: sympy/Quantum Computing.ipynb

max_issues_repo_path: sympy/Quantum Computing.ipynb

max_forks_repo_path: sympy/Quantum Computing.ipynb


text:
# 简单理解 SCF 中的 DIIS

> 创建时间：2019-10-23

这篇文档将会简单地叙述 GGA 为代表的 SCF DIIS。

DIIS 是一种 (几乎) 专门用于自洽场过程加速的算法。关于 DIIS 的算法与数学论述，这里并不作展开。这里推荐 C. David Sherrill 的笔记 [^Sherrill-note] 与 Psi4NumPy 的 Jupyter Notebook [^psi4numpy-note] 作为拓展阅读。

这篇笔记会借助 PySCF 的 DIIS 程序，对 Fock 矩阵进行外推。我们将描述在第 $t$ 步 DIIS 过程之后，如何更新第 $t+1$ 步的 Fock 矩阵。我们 **并不会** 从头写一个 DIIS 程序；这一方面是出于程序复杂性考虑，因为一般的 DIIS 程序应当要允许便利地增添、删除迭代过程中间的向量，能够处理解矩阵方程时会出现的线性依赖问题，并且要能保证一定的程序效率。另一方面，若对 DIIS 的更新过程有所了解，那么原则上我们已经理解了 DIIS 程序了，剩下的细节将只是时间与耐心上的问题。

```python
import numpy as np
import scipy
from pyscf import gto, dft, lib
from functools import partial

np.einsum = partial(np.einsum, optimize=["greedy", 1024 ** 3 * 2 / 8])
np.set_printoptions(5, linewidth=150, suppress=True)
```

## PySCF 的 DIIS 使用

### 分子体系定义与 DIIS 程序

首先，我们的分子体系是不对称的双氧水。

```python
mol = gto.Mole()
mol.atom = """
O  0.0  0.0  0.0
O  0.0  0.0  1.5
H  1.0  0.0  0.0
H  0.0  0.7  1.0
"""
mol.basis = "6-31G"
mol.verbose = 0
mol.build()

nao = nmo = mol.nao
nocc = mol.nelec[0]
nvir = nmo - nocc
so, sv, sa = slice(0, nocc), slice(nocc, nmo), slice(0, nmo)
```

为了简化程序，我们借助 PySCF 的 DFT 自洽场类 `scf_eng`，用来生成 Fock 矩阵和密度矩阵。

```python
scf_eng = dft.RKS(mol)
scf_eng.xc = "B3LYPg"
S = mol.intor("int1e_ovlp")
mo_occ = np.zeros(nmo)
mo_occ[:nocc] = 2
```

我们先简单地用下述程序展示 PySCF 中的 DIIS 是如何使用的。我们在后面会介绍 DIIS 类的具体的一些函数；这个用作演示用途的 DIIS 也时通过下述程序生成的。

下面的自洽场程序取自 pyxdh 文档 [^pyxdh-note]，但作了一些修改与简化。这个自洽场程序大体思路与 Szabo [^Szabo-Ostlund.Dover.1996] 第三章的叙述契合，也与 Psi4NumPy 和 PySCF 的不少演示程序比较相似。

与 Szabo 第三章叙述不太一样的地方有两处。其中一行代码是

```python
D = coef * D + (1 - coef) * D_old
```

这行代码仅仅是用来对 Naive SCF 作的修改。Szabo 的第三章可以称作是 Naive SCF，即单纯地将 Fock 矩阵对角化生成分子轨道，再得到密度代入 Fock 矩阵中。但这一行会将上一次迭代的密度 $D_{\mu \nu}^{t-1}$ 与这一次迭代的密度 $D_{\mu \nu}^{t}$ 混合，产生新的密度代入 Fock 矩阵中。这仅仅是为了防止 Naive SCF 振荡收敛，不会用于 DIIS 加速的算法。

另一行代码是

```python
F = func(diis=diis, F=F, C=C, mo_occ=mo_occ)
```

这行代码是用于指定 DIIS 的更新。方式在这份文档中，DIIS 的更新方式有

- `func_no_special`：Naive SCF，不引入 DIIS

- `diis_err_deviation`：通过迭代过程中的 Fock 矩阵差值 $\Delta F_{\mu \nu}^{t} = F_{\mu \nu}^{t} - F_{\mu \nu}^{t - 1}$ 更新 DIIS 状态

- `diis_err_gradient`：通过占据-非占 Fock 矩阵 $F_{ai}^{t}$ 更新 DIIS 状态

之所以用这种不太常见也不太直观的代码方式指定 DIIS 更新方式，单纯地是因为节省文档中的代码空间，避免代码冗余。

```python
def scf_process(func, coef=1.0, maxcycle=128):
    diis = lib.diis.DIIS()

    C = e = NotImplemented                                # Orbital (canonical) coefficient
    D = np.zeros((nao, nao))                              # Density in this iteration
    D_old = np.zeros((nao, nao)) + 1e-4                   # Density in last iteration
    count = 0                                             # Iteration count (1, 2, ...)

    while (not np.allclose(D, D_old)):                    # atol=1e-8, rtol=1e-5
        if count > maxcycle:
            raise ValueError("SCF not converged!")
        count += 1
        D_old = D
        F = scf_eng.get_fock(dm=D)                        # Generate initial Fock matrix from Density
        if count > 1:                                     # avoid the case: C = NotImplemented
            F = func(diis=diis, F=F, C=C, mo_occ=mo_occ)  # Different DIIS approaches
            # func_no_special    : nothing happens
            # diis_err_deviation : F = diis.update(F)
            # diis_err_gradient  : F = diis.update(F, scf_eng.get_grad(C, mo_occ))
        e, C = scipy.linalg.eigh(F, S)                    # Solve FC = SCe
        D = scf_eng.make_rdm1(mo_coeff=C, mo_occ=mo_occ)  # D = 2 * C(occ).T @ C(occ)
        D = coef * D + (1 - coef) * D_old                 # For convergence of original SCF
            # func_no_special: D = 0.3 * D + 0.7 * D_old
            # other cases    : nothing happens

    E_tot = scf_eng.energy_tot(dm=D)

    print("SCF Converged in     ", count, " loops")
    print("Total energy (B3LYP) ", E_tot, " a.u.")
```

### DIIS 方法加速效果

现在我们可以来看每种 DIIS 更新方式会产生的效果了。

**Naive SCF**

```python
def func_no_special(*args, **kwargs):
    return kwargs["F"]
scf_process(func_no_special, coef=0.3)
```

    SCF Converged in      51  loops
    Total energy (B3LYP)  -151.37754201912352  a.u.

**DIIS：Fock 矩阵差值**

```python
def diis_err_deviation(*args, **kwargs):
    diis, F = kwargs["diis"], kwargs["F"]
    return diis.update(F)
scf_process(diis_err_deviation)
```

    SCF Converged in      15  loops
    Total energy (B3LYP)  -151.37754323564036  a.u.

**DIIS：占据-非占 Fock 矩阵**

```python
def diis_err_gradient(*args, **kwargs):
    diis, F, C, mo_occ = kwargs["diis"], kwargs["F"], kwargs["C"], kwargs["mo_occ"]
    return diis.update(F, scf_eng.get_grad(C, mo_occ))
scf_process(diis_err_gradient)
```

    SCF Converged in      22  loops
    Total energy (B3LYP)  -151.37754323564045  a.u.

尽管从原理上，使用占据-非占 Fock 矩阵的方法应当更好；但在当前的体系下 Fock 矩阵差值的方法能更快地收敛。

我们能发现，若使用 PySCF 的 DIIS 类 (Psi4NumPy 的 [helper_HF.py](https://github.com/psi4/psi4numpy/blob/master/Self-Consistent-Field/helper_HF.py) 的 `DIIS_helper` 类也是相似的)，我们在实际进行 DIIS 时只需要相对于 Naive SCF 增加一行，非常方便地就可以加速收敛：

```python
F = diis.update(F)  # diis_err_deviation
```

或者

```python
F = diis.update(F, scf_eng.get_grad(C, mo_occ))  # diis_err_gradient
```

简单地说，这就是将 Fock 矩阵在每次迭代过程中，利用以前储存下来的 Fock 矩阵的信息进行再次更新。

## DIIS 细节

在这一段中，我们主要通过占据-非占 Fock 矩阵更新法的程序，对迭代 6 次时的 DIIS 状态进行较为细致的分析，并以此推导出第 6 次 DIIS 更新后的 Fock 矩阵。

第 6 次迭代时的 DIIS 类 `diis` 与更新前的 Fock 矩阵 `F_old` $F_{\mu \nu}^{t=6}$、更新后的 Fock 矩阵 `F` $\mathscr{F}_{\mu \nu}$ 的获得方式是通过限制迭代次数而给出的：

```python
diis = lib.diis.DIIS()

C = e = NotImplemented                                # Orbital (canonical) coefficient
D = np.zeros((nao, nao))                              # Density in this iteration
D_old = np.zeros((nao, nao)) + 1e-4                   # Density in last iteration
count = 0                                             # Iteration count (1, 2, ...)
F_old = NotImplemented                                # Variable in last iteration

while (not np.allclose(D, D_old)):                    # atol=1e-8, rtol=1e-5
    count += 1
    D_old = D
    F = scf_eng.get_fock(dm=D)                        # Generate initial Fock matrix from Density
    if count == 6:
        F_old = F.copy()
        F = diis.update(F, scf_eng.get_grad(C, mo_occ))
        break
    elif count > 1:                                   # avoid the case: C = NotImplemented
        F = diis.update(F, scf_eng.get_grad(C, mo_occ))

    e, C = scipy.linalg.eigh(F, S)                    # Solve FC = SCe
    D = scf_eng.make_rdm1(mo_coeff=C, mo_occ=mo_occ)  # D = 2 * C(occ).T @ C(occ)
```

这里补充一个程序细节，更新后的 Fock 矩阵 $\mathscr{F}_{\mu \nu}$ 也可以通过 `diis.extrapolate` 获得：

```python
np.allclose(F, diis.extrapolate().reshape(nmo, nmo))
```

    True

而 `diis.update` 除了给出更新后的 Fock 矩阵外，还将当前的 Fock 矩阵与误差信息更新入 `diis` 中，为下一次迭代的 DIIS 更新作准备。

### DIIS 储存内容

一般来说，DIIS 储存两部分内容：待外推信息 $p_I^t$ 与误差信息 $e_J^t$。

我们在 SCF 过程中使用 DIIS 的目的是借助以前若干步迭代过程中的信息，对 Fock 矩阵作外推，得到在当前迭代步 $t$ 下更好的 Fock 矩阵。因此，待外推信息 $p_I^t$ 是第 $t$ 次迭代过程计算得到的 Fock 矩阵 $F_{\mu \nu}^t$。

这里会有些诡异的地方是，待外推信息 $p_I^t$ 是单下标 $I$ 的向量，但原子轨道基组下的 Fock 矩阵 $F_{\mu \nu}^t$ 是双下标的矩阵。事实上，$p_I^t$ 在实践过程中就是将 $F_{\mu \nu}^t$ 压平成一维向量。待外推信息 $p_I^t$ 可以通过 `diis.get_vec` 给出，我们将这些向量整合为 `vecs` 变量中，其角标储存是 $(t, I)$：

```python
vecs = np.array([diis.get_vec(i) for i in range(diis.get_num_vec())])
```

我们记每次迭代的误差信息为 $e_J^t$。对于占据-非占 Fock 矩阵更新法而言，$e_J^t$ 即是分子轨道基组下的非占-占据 Fock 矩阵 $F_{ai}^t$。我们知道，对于 SCF 收敛之后的状态下，$F_{ai} = 0$；但这在自洽场迭代过程中，该量一般地并不是零，甚至说自洽场过程就是为了达成 $F_{ai} = 0$ 的目的也不为过。因此，$F_{ai}^t$ 的状况可以看作是自洽场是否收敛得较好的判标；于是我们定义 $e_J^t$ 为压平之后的 $F_{ai}^t$。

误差信息 $e_J^t$ 可以通过 `diis.get_err_vec` 给出，我们将这些向量整合为 `err_vecs` 变量中，其角标储存是 $(t, J)$：

```python
err_vecs = np.array([diis.get_err_vec(i) for i in range(diis.get_num_vec())])
```

我们指出，$p_I^t$ 与 $e_J^t$ 下标所指代的维度未必要是一样的。

```python
print(vecs.shape)
print(err_vecs.shape)
```

    (5, 484)
    (5, 117)

:::{note}

从上述的叙述与代码中，能看到我们只进行了 6 次迭代，其中迭代过程的误差信息与待外推信息只储存了 5 次 ($t$ 从 1 计数，外推信息的 $t \in [2, 6]$)。我们定义当前作为迭代次数的上标 $t$ 的集合是 $\mathscr{T} = \{2, 3, 4, 5, 6\}$；但在 PySCF 的 DIIS 类 `diis` 中，通过程序取出这些向量的指标时则应当使用 `0, 1, 2, 3, 4`。

我们知道，DIIS 为了进行外推，会储存许多待外推信息与误差信息；但对于大分子而言，这会占用许多内存空间。出于这个目的 (以及出于收敛性的考量)，DIIS 通常只会存比较少量地待外推信息与误差信息。PySCF 的 DIIS 一般只储存 6 次迭代过程的信息。这意味着，若我们进行了 15 次迭代，待外推矩阵至多也只会储存 6 个，其余的待外推信息或误差信息都会舍去。

为简化讨论，我们在这篇文档中不讨论如何舍去已经储存的待外推信息与误差信息。

:::

### DIIS 外推：理论

有了所有的待外推信息 $p_I^t$ 与误差信息 $e_J^t$ 后，我们可以作出外推结果 $\mathscr{p}_I = \mathscr{F}_{\mu \nu}$。上述公式中看似有问题的单下标转双下标可以通过互阵的 reshape 实现。

外推的含义是

$$
\mathscr{p}_I = \sum_{t \in \mathscr{T}} w_t p_I^t
$$

$\mathscr{T}$ 代表 DIIS 当前储存的每个外推信息对应的被迭代次数的集合，在这里恰好是从 2 开始的所有的被迭代次数。如果我们现在的迭代次数非常大，但只允许 DIIS 储存不多于 6 个待外推信息，那么求和指标 $t$ 的取值范围 $\mathscr{T}$ 将会舍去这些迭代次数，从而保持其集合元素数量 $|\mathscr{T}|$ 不超过 6。

我们人为地引入权重 $w_t$ 以归一条件：

$$
\sum_{t \in \mathscr{T}} w_t = 1
$$

如果我们假定待外推的信息 $p_I^t$ 与对应的误差信息 $e_J^t$ 呈线性关系，那么被外推的信息 $\mathscr{p}_I$ 的误差 $\mathscr{e}_J$ 应当满足

$$
\mathscr{e}_I = \sum_{t \in \mathscr{T}} w_t e_J^t
$$

我们希望误差 $\Vert \mathscr{e}_J \Vert_2^2$ 最小化，但同时又满足 $w_t$ 的归一化条件；那么我们通过 Lagrange 乘子法，构造以下损失函数

$$
\begin{align}
\mathscr{L} (\{w_t\}_{t \in \mathscr{T}}, \lambda) &= \Vert \mathscr{e}_J \Vert_2^2 + 2 \lambda \left( \sum_{t \in \mathscr{T}} w_t - 1 \right) \\
&= \sum_J \sum_{t \in \mathscr{T}} w_t e_J^t \cdot \sum_{s \in \mathscr{T}} w_s e_J^s + 2 \lambda \left( \sum_{t \in \mathscr{T}} w_t - 1 \right)
\end{align}
$$

我们现在定义

$$
B_{ts} = \sum_{J} e_J^t e_J^s
$$

那么损失函数可以写为

$$
\mathscr{L} (\{w_t\}_{t \in \mathscr{T}}, \lambda) = \sum_{t, s \in \mathscr{T}} w_t B_{ts} w_s + 2 \lambda \left( \sum_{t \in \mathscr{T}} w_t - 1 \right)
$$

对上述损失函数求取关于 $w_t$ 的偏导数，则得到

$$
\frac{\partial \mathscr{L}}{\partial w_t} = 2 \sum_{s \in \mathscr{T}} B_{ts} w_s + 2 \lambda
$$

我们显然是希望让损失函数对关于 $w_t$ 的偏导数为零；那么联立归一化条件 $\sum_{t \in \mathscr{T}} w_t = 1$，我们应当得到以下矩阵方程：

$$
\begin{align}
\begin{pmatrix}
  0 & 1 & 1 & \cdots \\
  1 & B_{t_0 t_0} & B_{t_0 t_1} & \cdots \\
  1 & B_{t_1 t_0} & B_{t_1 t_1} & \\
  \vdots & \vdots & & \ddots \\
\end{pmatrix}
\begin{pmatrix}
  \lambda \\ w_{t_0} \\ w_{t_1} \\ \vdots
\end{pmatrix}
=
\begin{pmatrix}
  1 \\ 0 \\ 0 \\ \vdots
\end{pmatrix}
\end{align}
$$

其中，$t_0, t_1, \cdots \in \mathscr{T}$ 是互不相同的指标。求解上述方程，就可以获得权重 $w_t$，进而给出 $\mathscr{F}_{\mu \nu} = \mathscr{p}_I = \sum_{t \in \mathscr{T}} w_t p_I^t$，达成目标。

### DIIS 外推：实现

首先，我们出，`diis` 的一个隐含变量 `diis._H` 储存的就是矩阵方程 LHS 的矩阵部分：

```python
A = diis._H[:diis.get_num_vec()+1, :diis.get_num_vec()+1]
A
```

    array([[ 0.     ,  1.     ,  1.     ,  1.     ,  1.     ,  1.     ],
           [ 1.     , 22.07074, -0.48488,  0.19048, -1.05937, -0.92243],
           [ 1.     , -0.48488, 21.30345,  0.23478,  1.63515, -0.41527],
           [ 1.     ,  0.19048,  0.23478,  1.35725, -0.18202,  0.38307],
           [ 1.     , -1.05937,  1.63515, -0.18202,  5.33326, -0.14979],
           [ 1.     , -0.92243, -0.41527,  0.38307, -0.14979,  1.18526]])

我们能很方便地构建上述矩阵的第 1 行以下、第 1 列以右的子矩阵 $B_{ts} = \sum_{J} e_J^t e_J^s$：

```python
np.einsum("tI, sI -> ts", err_vecs, err_vecs)
```

    array([[22.07074, -0.48488,  0.19048, -1.05937, -0.92243],
           [-0.48488, 21.30345,  0.23478,  1.63515, -0.41527],
           [ 0.19048,  0.23478,  1.35725, -0.18202,  0.38307],
           [-1.05937,  1.63515, -0.18202,  5.33326, -0.14979],
           [-0.92243, -0.41527,  0.38307, -0.14979,  1.18526]])

我们可以直接解上述的矩阵方程：

```python
b = np.zeros(diis.get_num_vec() + 1)
b[0] = 1
w = np.linalg.solve(A, b)
w = w[1:]
w
```

    array([0.05123, 0.02428, 0.31624, 0.13904, 0.46921])

那么我们可以通过 $\mathscr{F}_{\mu \nu} = \mathscr{p}_I = \sum_{t \in \mathscr{T}} w_t p_I^t$ 给出外推 Fock 矩阵 `F_ex`，并且与 `diis` 给出的外推的 Fock 矩阵 `F` 进行比较：

```python
F_ex = np.einsum("t, tI -> I", w, vecs).reshape(nmo, nmo)
np.allclose(F_ex, F)
```

    True

:::{tip}

在求解 DIIS 所给出的权重 $w_t$ 向量的过程中，会遇到线性依赖关系，或者说会遇到 $B_{ts}$ 数值上不满秩的情况。在这种情况下，求解矩阵方程可能会失败。

一种解决方案是，干脆关闭 DIIS，使用 Naive SCF 作最后的收尾工作。由于 DIIS 已经将电子态密度收敛到相当不错的状态了，因此应当能预期这种情况下 Naive SCF 可以正常地进行收敛。

另一种解决方式是对矩阵方程 $\mathbf{A} \boldsymbol{x} = \boldsymbol{b}$ 的矩阵 $\mathbf{A}$ 作对角化，并舍去其中绝对值极小的本征值与本征向量，求解一个子空间的线性方程组问题。这种解决方案应用在 PySCF 的 DIIS 程序中。

:::

### 对 Fock 矩阵差值方法的补充说明

Fock 矩阵差值方法的计算过程与占据-非占 Fock 矩阵方法的实现过程几乎是相同的。唯一的区别是：

- 占据-非占 Fock 矩阵方法 $e_J^t = F_{ai}^t$

- Fock 矩阵差值方法 $e_J^t = \Delta F_{\mu \nu}^{t} = F_{\mu \nu}^{t} - F_{\mu \nu}^{t - 1}$

[^Sherrill-note]: <http://vergil.chemistry.gatech.edu/notes/diis/diis.pdf>

[^psi4numpy-note]: <https://github.com/psi4/psi4numpy/blob/master/Tutorials/03_Hartree-Fock/3b_rhf-diis.ipynb>

[^pyxdh-note]: <https://py-xdh.readthedocs.io/zh_CN/latest/qcbasic/proj_xyg3.html>

[^Szabo-Ostlund.Dover.1996]: Szabo, A.; Ostlund, N. S. *Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Books on Chemistry)*; Dover Publications, 1996.


text:
# Notation: Converter Symbolic

```python
import mechkit
import numpy as np
import sympy as sp
import itertools
```

```python
np.set_printoptions(
    linewidth=140,
    precision=3,
    # suppress=False,
)
```

### Symbolic with numbers

```python
converter = mechkit.notation.ConverterSymbolic()
ones_tensor = np.ones((3, 3, 3, 3), dtype=sp.Symbol)
print(ones_tensor)
```

    [[[[1 1 1]
       [1 1 1]
       [1 1 1]]

      [[1 1 1]
       [1 1 1]
       [1 1 1]]

      [[1 1 1]
       [1 1 1]
       [1 1 1]]]

     [[[1 1 1]
       [1 1 1]
       [1 1 1]]

      [[1 1 1]
       [1 1 1]
       [1 1 1]]

      [[1 1 1]
       [1 1 1]
       [1 1 1]]]

     [[[1 1 1]
       [1 1 1]
       [1 1 1]]

      [[1 1 1]
       [1 1 1]
       [1 1 1]]

      [[1 1 1]
       [1 1 1]
       [1 1 1]]]]

```python
ones_mandel6 = converter.to_mandel6(ones_tensor)
print(ones_mandel6)
```

    [[1 1 1 sqrt(2) sqrt(2) sqrt(2)]
     [1 1 1 sqrt(2) sqrt(2) sqrt(2)]
     [1 1 1 sqrt(2) sqrt(2) sqrt(2)]
     [sqrt(2) sqrt(2) sqrt(2) 2 2 2]
     [sqrt(2) sqrt(2) sqrt(2) 2 2 2]
     [sqrt(2) sqrt(2) sqrt(2) 2 2 2]]

```python
ones_mandel9 = converter.to_mandel9(ones_tensor)
print(ones_mandel9)
```

    [[1 1 1 sqrt(2) sqrt(2) sqrt(2) 0 0 0]
     [1 1 1 sqrt(2) sqrt(2) sqrt(2) 0 0 0]
     [1 1 1 sqrt(2) sqrt(2) sqrt(2) 0 0 0]
     [sqrt(2) sqrt(2) sqrt(2) 2 2 2 0 0 0]
     [sqrt(2) sqrt(2) sqrt(2) 2 2 2 0 0 0]
     [sqrt(2) sqrt(2) sqrt(2) 2 2 2 0 0 0]
     [0 0 0 0 0 0 0 0 0]
     [0 0 0 0 0 0 0 0 0]
     [0 0 0 0 0 0 0 0 0]]

### Symbolic with letters

```python
def tensor(
    order=2, symbol="A", dim=3, latex_index=False, kwargs_symbol={}, indice_offset=0
):
    A = np.zeros((dim,) * order, dtype=sp.Symbol)
    for x in itertools.product(range(dim), repeat=order):
        index = "".join(map(str, map(lambda x: x + indice_offset, x)))
        if latex_index:
            index = "_{" + index + "}"
        A[x] = sp.Symbol(symbol + index, **kwargs_symbol)
    return A
```

```python
def make_it_hooke_symmetric(A, dim=3):
    for i in range(dim):
        for j in range(dim):
            for k in range(dim):
                for m in range(dim):
                    A[i, j, m, k] = A[i, j, k, m]
                    A[j, i, m, k] = A[i, j, k, m]
                    A[k, m, i, j] = A[i, j, k, m]
    return A
```

```python
def make_it_left_symmetric(A, dim=3):
    for i in range(dim):
        for j in range(dim):
            for k in range(dim):
                for m in range(dim):
                    A[j, i, k, m] = A[i, j, k, m]
    return A
```

```python
def make_it_right_symmetric(A, dim=3):
    for i in range(dim):
        for j in range(dim):
            for k in range(dim):
                for m in range(dim):
                    A[i, j, m, k] = A[i, j, k, m]
    return A
```

```python
def make_it_minor_symmetric(A, dim=3):
    tmp = make_it_left_symmetric(A)
    tmp = make_it_right_symmetric(A)
    return tmp
```

```python
tensor = make_it_minor_symmetric(tensor(order=4, indice_offset=1))
print(tensor)
```

    [[[[A1111 A1112 A1113]
       [A1112 A1122 A1123]
       [A1113 A1123 A1133]]

      [[A1211 A1212 A1213]
       [A1212 A1222 A1223]
       [A1213 A1223 A1233]]

      [[A1311 A1312 A1313]
       [A1312 A1322 A1323]
       [A1313 A1323 A1333]]]

     [[[A1211 A1212 A1213]
       [A1212 A1222 A1223]
       [A1213 A1223 A1233]]

      [[A2211 A2212 A2213]
       [A2212 A2222 A2223]
       [A2213 A2223 A2233]]

      [[A2311 A2312 A2313]
       [A2312 A2322 A2323]
       [A2313 A2323 A2333]]]

     [[[A1311 A1312 A1313]
       [A1312 A1322 A1323]
       [A1313 A1323 A1333]]

      [[A2311 A2312 A2313]
       [A2312 A2322 A2323]
       [A2313 A2323 A2333]]

      [[A3311 A3312 A3313]
       [A3312 A3322 A3323]
       [A3313 A3323 A3333]]]]

```python
tensor_mandel6 = converter.to_mandel6(tensor)
print(tensor_mandel6)
```

    [[A1111 A1122 A1133 sqrt(2)*A1123 sqrt(2)*A1113 sqrt(2)*A1112]
     [A2211 A2222 A2233 sqrt(2)*A2223 sqrt(2)*A2213 sqrt(2)*A2212]
     [A3311 A3322 A3333 sqrt(2)*A3323 sqrt(2)*A3313 sqrt(2)*A3312]
     [sqrt(2)*A2311 sqrt(2)*A2322 sqrt(2)*A2333 2*A2323 2*A2313 2*A2312]
     [sqrt(2)*A1311 sqrt(2)*A1322 sqrt(2)*A1333 2*A1323 2*A1313 2*A1312]
     [sqrt(2)*A1211 sqrt(2)*A1222 sqrt(2)*A1233 2*A1223 2*A1213 2*A1212]]

```python
tensor_mandel9 = converter.to_mandel9(tensor)
print(tensor_mandel9)
```

    [[A1111 A1122 A1133 sqrt(2)*A1123 sqrt(2)*A1113 sqrt(2)*A1112 0 0 0]
     [A2211 A2222 A2233 sqrt(2)*A2223 sqrt(2)*A2213 sqrt(2)*A2212 0 0 0]
     [A3311 A3322 A3333 sqrt(2)*A3323 sqrt(2)*A3313 sqrt(2)*A3312 0 0 0]
     [sqrt(2)*A2311 sqrt(2)*A2322 sqrt(2)*A2333 2*A2323 2*A2313 2*A2312 0 0 0]
     [sqrt(2)*A1311 sqrt(2)*A1322 sqrt(2)*A1333 2*A1323 2*A1313 2*A1312 0 0 0]
     [sqrt(2)*A1211 sqrt(2)*A1222 sqrt(2)*A1233 2*A1223 2*A1213 2*A1212 0 0 0]
     [0 0 0 0 0 0 0 0 0]
     [0 0 0 0 0 0 0 0 0]
     [0 0 0 0 0 0 0 0 0]]

max_stars_repo_name: jonashund/mechkit

max_issues_repo_name: jonashund/mechkit

max_forks_repo_name: jonashund/mechkit