text:
### 1.5 Decision Theory

##### 이 장의 목적

   * 1.2 절에서 확률이론을 사용하여 불확실성에 대해서 어떻게 처리 하는지 수학적 기법을 소개 하였다.
   * 결정 이론과 확률론이 결합하여, 불확실성이 포함된 상황에서 최적의 결정을 할 수 있으며, 패턴인식에서 그러한 상황을 만난다.
   * 따라서 Decision Theory 에 확률이론이 어떻게 적용하는 기법들에 관해서 논의한다.

##### 주어진 과제

   * 입력벡터 ${\bf x}$ 와 출력(목표) 벡터 ${\bf t}$ 를 보자.
   * 과제는 새로운 값 ${\bf x}$ 가 주어질때, ${\bf t}$ 를 에측 할 수 있는지???

###### 상화에 따른 과제의 성격

   * Regression 문제
        * ${\bf t}$ 는 연속변수에 대한 예측 문제

   * 분류문제
       * ${\bf t}$ 는 라벨이 될 것이다.

   * 결합확률분포 $p({\bf x}, {\bf t})$ 는ㄴ 이러한 변수들의 불확정성에 대해서 기술을 한다.
   * $p({\bf x}, {\bf t})$ 를 결정하는 것은 학습 데이타로 부터 추론 하는 것이며, 일반적으로 이책의 많은 부분에서 해법을 제시하는것은 어렵다.

   * 그러나 실제적인 문제에서, ${\bf t}$ 에 대해서 구체적인 예측을 할 수 있어야만 하고, 더욱 일반적으로는 ${\bf t}$ 가 취해야 할 값에 대해 이해한 바탕으로 구체적인 행동을 하여야 한다.
   * 이러면이 decision theory 의 주제이다.

##### 메디컬 진단 예

   * X-ray 사진이 있으며, 사진으로 환자가 암인지 아닌지 판정해야 한다.

   * 입력 벡터 ${\bf x}$ 는 이미지에서 픽셀의 강도가 되며, 출력 ${\bf t}$ 는 암의 존재 여부가 된다.
   * 이제 분류를 $\mathcal C_1$ 을 암이 있는 이미지, $\mathcal C_2$ 를 암이 없는 이미지 라고 하자.
   * ${\bf t}$ 는 이진값으로 $t = 0$ 는 ${\mathcal C_1}$ 으로 $t = 1$ 는 ${\mathcal C_2}$ 로 대응 할 수 있다.

   * 이런식으로 라벨링 하는 방법은 확률모델에서 매우 편리한 방법인지는 나중에 보여 줄것이다.

##### Decision Theory 의 수학적 형태.

   * 메디컬 진단의 예에서와 같이, 일반적인 추론에 대한 문제는 결합확률분포 $p({\bf x}, {\mathcal C}_k)$ 이나 $p({\bf x}, {\bf t}) 를 결정하는 문제이다.
   * 이런 수학적 형식이 유용하고 유익한 양이 될지라도, 결국에는 환자를 치료해야 할지 말지를 결정해야 하며, 이 결정이 적절한 선택이 되기를 희망한다.
   * 이것인 decision step 이며, 주어진 환경에서 최적의 결정을 하는지에 대한 decision theory 의 주제이다.
   * 일단 추론문제를 풀고나면 decision stage 는 간단하고 매우 평이하다.

##### decision theory 의 key idea

   * decision theory 의 이책의 나머지에서 요구하는 주요 아이디어를 소개.
   * X-ray 이미지에서 이미지에 연결된 두가지 분류중에 하나를 결정해야 한다.
   * 베이지안 이론을 적용한다.
$$p(C_k \mid {\bf x})=\dfrac{p({\bf x}\mid C_k)p(C_k)}{p({\bf x})} \qquad{(1.77)}$$

   * $p(C_k)$ 를 $C_k$ 에 대한 사전확률로 해석하고, $p(C_k \mid {\bf x})$ 를 관련있는 사후확률로 해석한다.
   * 그러므로, $p(C_1)$ 을 암이 있는 환자의 확률이고, $p(C_1 \mid {\bf x})$ 를 X-ray 가 주어진 상황에서 암이 있는 확률이다.

   * 목적은 잘못된 분류를 최소화 하는 것이므로, 직관적으로 사후확률이 높은 분류를 갖는 것으로 선택할 것이다.
   * 이런 직관이 맞는것이지, 그리고 decision 을 하기 위한 일반적인 조건이 무엇인지 대한 논의를 할 것이다.

#### 1.5.1 오분류 최소화 (Minimizing the misclassification rate)

##### Decision Region

   * 각 변수 ${\bf x}$ 에 대해서 가능한 분류 하나씩 대응한다.
   * 이것은 입력 변수를 각각 영역 $\mathcal R_k$ 에 대응시킨다라고 하고, $\mathcal R_k$ 를 Decision Region 이라고 한다.
   * $\mathcal R_k$ 에 있는 모든 점은 분류 $\mathcal C_k$ 에 대응된다.

   * Decision Region 사이에 있는 경계를 decision boundaries 또는 decision surface 라고 한다.
   * 각 decision 영역은 인계될 필요는 없지만, 어떤 분리된 영역을 포함 할 수 있다.
   * 이장의 후에, decision region 과 decision boundary 의 예를 만나게 될것이다.

##### 잘못분류  된 경우 확률

   * 최적의 결정을 하기 위해서는 암문제 처럼, 모든 경우를 먼저 고려 하여야 한다.
   * 잘못된 결정은 $C_1$ 에 속한 벡터를 $C_2$ 로 분류하거나, 또는 그 반대로 할 경우에 발생한다.
   * 이것이 발생활 확률은 다음과 같다.

 $$p(mistake)=p(x\in \mathcal R_1, C_2) + p(x\in \mathcal R_2, C_1)=\int_{\mathcal R_1}p({\bf x}, C_2)d{\bf x} + \int_{\mathcal R_2}p({\bf x}, C_1)d{\bf x} \qquad{(1.78)}$$

#### 각 $x$ 를 두 클래스 중에 연결하는 방법

   * 각 $x$ 를 두 클래스 중에 하나에 부여는 자유롭게 할 수 있다. 분명하게 $p(mistake)$ 를 최소화 하기 위해서는 각 $x$ 는 식 (1.78) 에서 작은 값을 갖는 클래스에 부여 하여야 한다.
   * 따라서 주어진 $x$ 에 대해서 $p({\bf x}, \mathcal C_1) \gt p({\bf x}, C_2)$ 이면 $x$ 는 class $C_1$ 에 부여 해야 한다.

   * 확률곱 에 의하여 $p({\bf x}, C_k) = p(C_k \mid {\bf x})\, p({\bf x})$ 가 된다.

   * $p({\bf x}, \mathcal C_1) \gt p({\bf x}, C_2)$ 는 $p(C_1 \mid {\bf x})\, p({\bf x}) \gt p(C_2 \mid {\bf x})\, p({\bf x})$ 가 되고, 양쪽 항에 공통이 $p({\bf x})$ 를 제거하면

   * 결국 mistake 를 최소화 하는 것은 각 ${\bf x}$ 를 사후 확률 $p(C_k \mid {\bf x}) 가 큰쪽 클래스로 부여하는 것이다.

   * 이런 결과,  단일 변수 x, 두개 클래스에 대해서 그림 1.24 에 표시 되어 있다.

   * 그림 1.24 x 에 대한 두 클래스 각각에 대한 결합 확률 $p(x, C_k)$ 그래프.

   * 결정경계 $x = {\hat x}$ 를 포함 $x$ 에 대해 구분된 두 클래스 각각에 대한 결합확률 $p(x, C_k)$ 도식화.
   * $x \geq {\hat x}$ 는 $C_2$ 로 분류하고, decision region 은 $\mathcal R_2$ 에 속하게 된다.
   * $x \lt {\hat x}$ 는 $C_1$ 로 분류하고, decision region 은 $\mathcal R_1$ 에 속하게 된다.
   * 에러는 불루, 그린, 레드 영역에서 발생한다. 즉 $x \lt {\hat x}$ 에서 클래스 $C_2$ 가 $C_1$ 로 잘못 분류된다.(레드와 그린이 합친영역)
   * 반대로, $x \geq {\hat x}$  에서 클래스 $C_1$ 을 $C_2$ 로 잘못 분류된다(불루 영역)
   * decision boundary ${\hat x}$ 가 이동하면, 불루하고 그린이 합쳐진 영역은 변하지 않으나($C_1$ 을 $C_2$ 로 오인), 레드 영역은 변하게 된다.
   * ${\hat x}$ 의 최적의 선택은 곡선 $p(x, C_1)$ 과 $p(x, C_2)$ 가 교차하는 점이고, 이 영역에선  ${\hat x} = x_0$ 이고 레드 영역이 사라진다.
   * 이 점에서 misclassification 비율이 최소가 되는 것과 같다. 그리고 $x$ 의 각 값에 대해서 높은 사후확률 $p(C_k \mid x)$ 를 갖는다.

#### 클래스를 K 개로 확장하면

   * K 개 클래스의 일반적인 경우, correct 되는 확률을 최대화 하는것은 다음과 같다.
  $$p(correct)=\sum_{k=1}^{K}p(x \in \mathcal R_k, \mathcal C_k) = \sum_{k=1}^{K}\int_{\mathcal R_k}p(x, \mathcal C_k) d{\bf x} \qquad{(1.79)}$$

   * 각 ${\bf x}$ 가 $p({\bf x}, C_k)$ 가 최대가 되는 클래스에 부여되는 그런 영역 $R_k$ 에 선택되면, 식 1.79 는 최대가 된다.
   * $p({\bf x}, C_k) = p(C_k \mid {\bf x})\,p({\bf x})$ 가 되고, $p({\bf x})$ 가 공통 이면, ${\bf x}$ 는 사후확률 $p(C_k \mid {\bf x})$ 가 최대가 되는 클래스에 부여되어야 한다.

#### 1.5.2 Minimizing the expected loss (기대손실 최소화)

##### 의료 진단 문제고 다시 가서

   * 많은 경우 단순히 오분류 수를 최소화 하는 것 보다는 더 복잡하다.
   * 의료 진단 문제로 다시 가서
   * 만일 만일 암이 아닌 환자가 잘못 진단하여 암으로 진단하면(첫번경우), 결과적으로 환자의 고통과 추가 조사의 필요성이 있을 것이다.
   * 반대로, 암환자가 건강한것으로 진단을 내면 결국은 치료부족으로 조기 사망이 될 수도 있다(두번째 경우).
   * 그러므로, 두 가지 경우 잘못진단 결과는 크게 다를 수 있다.
   * 두번째 경우가 첫번째 경우보다 비용이 더 들어도, 첫번쨰 경우보다 실수를 더 줄이는 것이 좋다.

#### loss function

   * 이런 문제는 $loss\,function$ 도입함으로써 형식화 할 수 있다.
   * $loss\,function$ 는 종종 $cost\,function$ 이라고 한다.
   * 이 함수는 가능한 decision 이나 행위를 할때 발생할 손실의 전반적이고 유일한 척도 이다.
   * 목적은 발생할 전체적인 손실을 최소화 하는 것입닏다.
   * 일부는 $loss\,function$ 대신에 $utility\,function$ 을 도입하여 이값이 최대가 되도록 한다.
   * $utility$ 가 단순히 $loss$ 를 음수화 한다면, 개념적으로 같다.

##### loss function 형식화

   * 새로운 값 ${\bf x}$ 가 있고, 그것의 진짜 클래스는 $C_k$ 이다.
   * 이 ${\bf x}$ 를 클래스 $C_j$ 에 부여했다. ($j$ 는 반드시 $k$ 와 같을 필요가 없다.)
   * 이럴경우 어떤 수준의 손실이 발생하는데 $L_{kj}$ 로 표시 하도록 한다.
   * $L_{kj}$ 를 $loss\,matrix$ (손실행렬) 에서 $k,j$ 요소로 나타난다.
   * 암진단 예에서 보면 다음 그림 1-25 형식이 된다.

 그림 1-25

    * row : true class
    * column : decision 기준에 의하여 부여된 class

   * 그림 1-25 에서 보면 correct decirion 이 발생하면 손실은 없다. (암환자는 암 진다, 건강한 사람은 건강한것으로 진단)
   * 건강한사람을 암진단하면 1, 암환자를 건강한 사람으로 진단 1000.

##### loss function 최적화 해법

   * $loss\,function$ 을 최소화 하는것이 최적해 이다.
   * 입력 벡터 ${\bf x}$ 가 주어질때고, 진짜 클래스 의 불확정 정도는 결합확률 $p({\bf x}, \mathcal C_k)$ 로 표현한다.
   * 평균 손실을 최소화 하길 원한다.
   * 평균은 다음과 같이 계산한다.

 $$\mathbb E[L] = \sum_{k}\sum_{j}\int_{R_j}L_{kj}\,p({\bf x}, C_k)\,d{\bf x} \qquad{(1.80)}$$

   * 각 ${\bf x}$ 는 독립적으로 decision 영역 $\mathcal R_j$ 에 부여된다.
   * 목적은 식 손실 기대값 1.80 이 최소화 하도록 $R$ 을 선택하는 것이다.
   * 식 1.80 에서, ${\bf x}$ 에 대해서, $\sum_{k}L_{kj}\,p({\bf x}, C_k)$ 를 최소화 해야 한다.

   * $p({\bf x}, C_k) = p(C_k \mid {\bf x})\,p({\bf x})$ 를 이용하고, $p({\bf x})$ 공통을 제거한다.

   * 기대손실을 최소화 하는 decision rule 은 각 ${\bf x}$ 에 대해서 클래스 $j$ 에 부여되는 손실 양은

$$\sum_{k} L_{kj} \,p(C_k \mid {\bf x}) \qquad({1.81})$$

   * 식 1.81 이 최소화 하는 것이어야 한다.
   * 이것은 잘 알다시피, 클래스에 대한 사후 확률을 알면 명백하게 평범한것이다.

#### 1.5.3 The Reject Option (리젝트 옵션)

   * 분류에러는 사후확율 $p(C_k \mid {\bf x})$ 가 1 보다 현저하게 작거나, 또는 결합분포 $p({\bf x}, C_k)$ 와 유사한 값을 갖는 , 입력 space 영역 에서 발생한다.
   * 이런 영역은 클래스에 속할지 불확실한 경우가 있는 영역이다.
   * 일부 응용시스템에서는 분류 결정을 해야 하는 경우 낮은 에러를 내기 어려운 경우에 decision 을 하지 않는다. 이러 것을 $rejection\,option$ 이라고 한단.
   * 예를 들어, 의료예제에서, X-ray 이미지에서 의심의 여지가 없는 것은 자동 분류를 하고, 애매한 경우에 전문가에게 넘긴다.
   * 이런 경우 임계값 $\theta$ 를 도입하여, 사후확률 $p(C_k \mid {\bf x})$ 의 최대값이 $\theta$ 보다 작거나 같은 그런 ${\bf x}$ 에 대해서 reject 할 수 있도록 한다.

   * 그림 1-26에서 단일 연속변수 $x$ , 2가지 클래스에 대해서 설명한다.

   * 그림 1-26 reject region

   * $x$ 에 대해서 사후 확률 이 임계값 $\theta$ 보다 같거나 작은 경우 reject 한다.

   * $\theta = 1$ 이면 모든 경우가 reject 한다.
   * $\theta \lt \frac{1}{K}$ 가 있어서 reject 되는 것이 없는 클래스를 K 로 잡는것으로 하여, reject 를 $\theta$ 에 의해서 조정하게 된다.
   * 이것을 reject decision 이 발생할때, 발생할 손실을 반영하여, loss matrix 를 주어지고, expected loss 를 최소화 하는 기준으로 reject 기준을 확장시켜 나갈수 있다.

### 1.5.4 Inference and decision (추론 과 결정)

   * 분류문제를 두 단계롤 나눈다. 추론단계 (inference stage) 와 decision stage(결정단계)
   * 추론 단계에서는 학습데이터로 $p(C_k \mid {\bf x})$ 를 결정
   * 계속해서 decision stage(결정단계) 에서 사후확률을 이용하여 최적의 클래스를 부여한다.

   * 또다른 접근 방법은 두개의 문제를 한꺼번에 해결하고, 간단하게 입력 ${\bf x}$ 를 decision 에 직접 맵핑 하는 함수를 학습하는 것이다. 이런 함수를 $discriminant\,function$ (판별함수) 이라고 한다.

##### 실제문제에서 decision problems 해결에서 접근방식

   * decision problems 해결에서 3가지 접근 방식을 정의 할 수 있다. 3가지 다 실제 문제 해결에서 사용한다. 복잡도가 큰것 부터 설명한다.
   * 문제는 주어진 ${\bf x}$ 가 어느 클래스 $C_k$ 에 속하는지 추론하는 방법을 푸는데는 3가지의 접근 방법을 사용

##### generative models

   * 먼저 각 클래스 $C_k$ 에 대한 클래스 조건 밀도를 $p(C_k \mid {\bf x})$ 와 사전 확률 $p(C_k)$ 를 구하고 나서 이것을 이용해서 사후 확률 $p(C_k \mid {\bf x})$ 를 구한다.
   * 베이즈 이론 을 사용하요 사후 확률을 구한다.
   $$p(C_k|{\bf x})=\dfrac{p({\bf x}|C_k)p(C_k)}{p({\bf x})} \qquad{(1.82)}$$
   * 분모 는 다음과 같다.
   $$p({\bf x}) = \sum_{k} p({\bf x}) \mid C_k)\, p(C_k) \qquad({1.83})$$

   * 직접 결합확률을 사용하여 구하는 경우도 있다 결합확률은 $p({\bf x}, C_k)$ 이다.

   * 사후 확률을 찿은후에는 decision theory 를 이용하여 각 신규 ${\bf x}$ 마다 해당되는 클래스를 결정합니다.
   * 이런 방법의 입력 에 대한 분포를 가지고 있어서, 그 분포로 부터 샘플링이 가능하기 때문에 generative model 이라고 한다.

#### discriminativie models

   * 먼저 사후 확률을 $p(C_k \mid {\bf x})$ 를 직접 구하고 나서, decision theory 를 이용하여 신규 ${\bf x}$ 에 대해서 클래스를 찿는다.
   * 이런 접근은 사후 확률로서 직접 찿는 모델이라서 discriminative models 라고 한다.

#### discriminant function

   * discriminant function 이라고 불리는 $f(x)$ 를 구한다. 이 함수는 각 입력 $x{\bf x}$ 에 대해서, 직접 클래스에 매핑 하는 함수 이다.
   * 예를 들면, 2 클래스 문제에서, $f(.)$ 는 binary 함수가 될 것이고 $f = 0$ 는 $C_1$ 에, $f = 1$ = $C_2$ 에 직접 설정한다.
   * 이 경우에는 확률을 사용하지 않는다.

####  상대적인 장점들..

###### Generative models

   * ${\bf x}$ 와 $C_k$ 의 결합 확률이므로, 많이 요구한다.
   * 많은 시스템에서, ${\bf x}$ 는 고차원이므로, $p({\bf x} \mid C_k)$ 를 합리적으로 정학하게 구하기 위해서 많은 학습데이터를 요구한다.
   * 사전확률 $p(C_K)$ 는 종종 각 클래스의 데이터 포인트의 비로서 추정하기도 한다.
   * 이 접근 방법의 장점은 식(1.83) 데이터 $p({\bf x})$ 의 marginal density 를 구할 수 있다는 것이다.
   * 이 방법은 낮은 확율을 갖고, 낮은 정확도의 예측을 하는 신규데이터를 찿아내는데 유용하다. 이것은 outlier detection or nvelty detection 으로 알려져 있다.
   * 그러나 단지 분류 결정만 하기에는, 너무 많은 데이터와, 컴퓨터 자원을 사용한다.
   * 사후확률 $p(C_k \mid {\bf x})$ 를 구하려고 하는데, 결합확률 $p({\bf, x}, C_k)$ 를 구하려고 많은 자원을 사용한다.

##### Discriminative models

   * 클래스 조건 확률밀도는 많은 것을 포함하고 있지만, 사후확률에는 별로 영향이 없는 것도 있다.
   * 그림 1-27 에 이런 내용을 표시한다.

   * 그림 1-27 클래스 조건 확률과 입력조건 확률 비교

   * 단 변수 $x$ 를 갖는 두 클래스 의 class 조건 확률밀도(왼쪽) 과 대응되는 사후확률밀도 (오른쪽)
   * 좌변의 $p(x \mid C_1)$ 은 푸른색 그래프이고, 우변의 사후확률에 영향을 주지 못한다... 왜??/
   * 우변의 녹색 수직선은 $x$ 의 결정경계를 나타내며 거기서 최소 오분류율을 나타낸다.

   * generative models 과 discriminative models 의 기계학습에서 상대적 장점들을 찿는것에 그리고 그것들을 합치는 방법에 대해서 많은 흥미가 있다.

##### Discriminative function

   * 아주 간단한 방법은 학습데이터로 부터 discriminative function 을 직접구하는것이다.
   * 이 함수는 각 ${\bf x}$ 에 직접 클래스를 연결한다. 따라서 추론과 결정단계가 하나의 학습문제로 되버린다.
   * 그림 1.27 에서 x 가 수직 녹색선을 찿는것과 같다. 즉 오분류의 최소 확률를 주는 결정경계이기 때문이다.
   * 그러나, 이 접근 방법은 더이사 사후 확률을 사용하지 않는데, 사후 확률을 사용해야만 하는 강력한 이유들이 있다.

#### 사후 확률을 사용해야만 하는 이유

###### 리스크 최소화 (Minimizing Risk)

   * loss matrix 가 시간마다 변화는 그런 시스템을 고려하자 (financial system 에서 나타난다)
   * 사후 확률을 안다면 1.18을 이용하여 적절하게 risk decision 을 최소화 하는 기준을 시시각각 갱신 할 수 있다.
   * discriminative funtion 만 있으면, loss matrix 는 다시 학습데이타로 가서 분류문제를 풀려고 할것이다.

###### Reject Option

   * 사후 확률이 오분류를 최소화 하는  rejection 기준을 결정하도록 한다. 또한 더해서 일반적으로 rejected data point 에 대한 기대손실을 알려준다.

###### Compensating for class priors (클래스들의 사전 분포에 대한 보상)

   * 다시 X-ray 문제로 돌아가서
   * 자동 스크린닝 시스템을 의해서, 학습데이터로 사용하기위해서 일반인으로 부터 많은수의 X-ray 이미지를 얻었다고 하자.
   * 암 경우는 일반인에 비해서 많지 않으므로, 그러한 데이타를 학습으로 사용한다면, 암 클래스에 너무 적은 빈도 때문에 어려움에 닥칠것이다.
   * 예를 들어, 모든 점을 정상 클래스에 부여하는것은 이미 99.9%의 정확도를 가지고 있으며, 이런 상식적인 해를 피하기는 어렵다.
   * 또한 많은 데이타 셋트가 아주 적은 수의 암 관련 이미지를 가지고 있기 때문에 학습 알고리즘은 암을 가지고 있는 이미지에 잘 노출이 되지 않아서, 일반화가 잘 되지 않는다.
   * 각 클래스에 같은수의 균형잡힌 데이터 셋트가 더 정확한 모델을 구현 할 수 있다.
   * 이 경우에는 수정된 효과를 학습데이터에 반영하여야만 한다.
   * 수정된 학습데이터를 사용하고, 사후확률에 대한 모델을 알고 있다면,
   * 베이즈 이론 1.82 에 의해서, 사후확률은 사전확률에 비례하고 있으며, 각 클래스에 대한 데이터 점들의 비율 로서 해석 될 수 있다.
   * 따라서, 인위적인 균형잡힌 데이터로 부터 얻은 사후확률을 얻을수 있으므로 (단지 카운팅으로 가능),
   * 먼저 데이터를 클래스 비율로 나눈다음
   * 모델에 적용하기 위한 전체수에 클래스 비를 곱해준다.
   * 마지막으로 전체 사후 확율 합이 1이 되도록 정규화한다.
   * 이과정은 사후확률 성질을 이용하는 대신에 직접 discriminative function 을 훈련한다면 적용할수 없다.

###### Combining models

   * 복잡한 시스템에서, 각 분리된 모듈에서 담당하도록 문제를 여러 하부 문제로 나누기를 원한다.
   * 예를 들면, 의료진단 문제에서, X-ray 이미지 이외에, 혈액 테스트 정보를 가지고 있다.
   * 그럼 이질적인 정보를 모두 모아 한꺼번에 하기 보다는 각각 X-ray 이미지 를 담당하고, 다른것은 혈액정보를 담당하게 함으로써 효과적으로 모델을 구성 할 수 있다.
   * 2 모델의 각각이 클래스별로 사후확률을 주면, 확률의 법칙을 사용하여 확률을 합칠수 있다.
   * 간단한 방법은 다음과 같다.

###### 확률을 합치는 방법

   * 입력자료를 X-ray 이미지는 ${\bf x}_I$, 혈액정보는 ${\bf x}_B$ 는 독립적이다고 하면

$$p({\bf x}_I, {\bf x}_B \mid C_k) = p({\bf x}_I \mid C_k)\,p({\bf x}_B \mid C_k) \qquad({1.84})$$

  * 사후확률은 다음과 같이 된다.

  $$\begin{align} p(C_k \mid {\bf x}_I, {\bf x}_B) \propto p({\bf x}_I, {\bf x}_B \mid C_k)\,p(C_k) \\\\ \propto p({\bf x}_I \mid C_k)\,p({\bf x}_B \mid C_k)\,p(C_k) \\\\ \propto \dfrac{p(C_k \mid {\bf x}_I)\,p(C_k \mid {\bf x}_B)}{p(C_k)} \end{align} \qquad{(1.85)}$$

   * 사전확률 $p(C_k)$ 가 필요하며, 클래스에서 데이터 포인트의 비율로 부터 추정할 수 있다. 그런다음, 전체합이 1이 되도록 정규화 한다.
   * 1.84 가정은 naive Bayes model 에서 적용 할 것이다.
   * 결합확률 $p({\bf x}_I, {\bf x}_B)$ 는 일반적으로 이 모델에서는 factorize 되지 않는다.
   * 이것은 마지막장에서 조건부 독립가정 1.84를 요구하지 않는 데이타 결합하는 모델에 대해서 논의 할 것이다.

### 1.5.5 Loss function for regression

   * 분류 문제에 대해서 decision theory 를 논의 하였다.
   * 다시 앞에서 논의하였던 커브피팅에 대해서 regression problems 로 다시 논의 한다.
   * decision stage 는 각 입력 변수 ${\bf x}$ 에서 값 $t$ 의 추정함수 $y({\bf x})$를 선택하는 것을 포함한다.

   * 손실함수 $L(t, y({\bf x}))$ 라고 하면
   * 손실 평균 또는 기대값은 다음과 같이 한다.

   $$\mathbb E[L] = \iint L(t, y({\bf x}))p({\bf x}, t)\,d{\bf x}dt \qquad{(1.86)}$$

   * 앞에서 논의한 sum-of-square 대신에 loss function 을 사용하였다.
   * sum-of-square 경우에는 최소가 되는 w 를 구하였는데, 여기서는 그냥 함수 y를 변수로 갖는 범함수(functional) 로 정의 하였다.
   * 따라서 $L$ 은 $y$ 의 함수이므로 $L$ 이 최소가 되는 함수 $y$ 의 형태를 묻게 된다.

   * 회귀문제에서 손실함수의 일반적은 것은 다음과 같이 제곱 손실형태로 주어진다.
   * $L(t, y({\bf x})) = \left \{ y({\bf x}) - t \right \}^2$

##### 기대손실의 최소화

   * 기대 손실은 손실함수가 주어진 형태에서 다음과 같이 된다.

   $$ \mathbb E[L] = \iint \left \{y({\bf x})-t \right \}^2 p({\bf x}, t)\,d{\bf x}dt \qquad{(1.87)}$$

  * 우리목적은 $\mathbb E[L]$ 을 최소화 하는 함수 $y$ 를 선택하는것이다. $\mathbb E[L]$ 은 범함수 (functional) 이다.

  * 변분원리를 적용하여

$$\dfrac{\delta \mathbb E[L]}{\delta y({\bf x})}=2\,\int\left \{y({\bf x})-t\right \}\,p({\bf x}, t)\,dt=0 \qquad{(1.88)}$$

   * $y({\bf x})$ 에 대해서 풀어보면

$$y({\bf x})= \dfrac{\int t\,p({\bf x}, t)\,dt}{p({\bf x})}=\int t\,p(t \mid {\bf x})dt=\mathbb E_t[t \mid {\bf x}] \qquad{(1.89)}$$

   * 이것은 ${\bf x}$ 가 주어진 에서 $t$ 의 평균값을 구하는 것이고, regression function 으로 알려져 있다.
   * 이것은 그림 1.28 에 표시되어 있다.

   * 그림 1-28 regression function $y(x)$

   * 조건 확률 분포 $p(t \mid x)$ 가 주어진 상태에서, expected squared loss 를 최소화 한다.

##### 다변량 타켓 변수 인 경우

   * $t$ 를 다변량으로 확장하면, t 는 벡터가 된다. ${\bf t} $
   * 이 경우 최적해는 ${\bf y}({\bf x}) = \mathbb E_t[{\bf t} \mid {\bf x}]$ 가 된다.

##### 최적해를 구하는 다른방법

   * 위 결과를 조금 다른 방법으로 유도함. 회귀문제에 대한 기본성질을 밝힐것이다.
   * 조건 기대치에 대한 최적해에 대한 풀이 방법을 생각해서,
   * 1.87 의 적분안에 있는 제곱항을 다음과 같이 전개한다.

 $$\left\{y({\bf x})-t\right\}^2 = \left\{y({\bf x})-\mathbb E[t\mid{\bf x}]+\mathbb E[t\mid{\bf x}]-t\right\}^2\\ = \left \{y({\bf x})-\mathbb E[t\mid {\bf x}] \right \}^2 + 2\left \{y({\bf x})-\mathbb E[t\mid {\bf x}] \right \}\left\{\mathbb E[t\mid {\bf x}]-t\right \}+ \left \{\mathbb E[t\mid {\bf x}]-t \right \}^2$$

   * 표기를 깔끔하기 하기 위해서, $\mathbb E_t[t \mid {\bf x}]$ 를 $\mathbb E[t \mid {\bf x}]$ 로 표시 할 것이다. (1.89 참조)
   * $y({\bf x})$ 는 예측할려고 하는 값이며, 모델에 의해서 결정된다.
   * 식 1.89 를 참조하면, $\mathbb E[t \mid {\bf x}]$ 는 t 의 평균값이다.

   * 이 식을 1.87 에서 loss function 에 넣고 t 로 적분하면 cross-term 은 없어지고 ( 없어지게 적당히 $\mathbb E[t \mid {\bf x}]$ 를 잡는다), 다음과 같은 형태로 된다.

$$\mathbb E[L] = \int \left\{y({\bf x})- \mathbb E[t\mid{\bf x}]\right\}^2 p({\bf x})\, d{\bf x} + \int \left\{\mathbb E[t \mid {\bf x}]- t \right\}^2 p({\bf x})\,d{\bf x} \qquad{(1.90)}$$

   * 식 1.90 의 두번째 항은 $\int var[t\mid {\bf x}]\,p({\bf x})\,d{\bf x} $ 과 같다.

   * 첫번째 항 에서 $y({\bf x})$ 는 우리가 찿을려고 하는 식이며, $\mathbb E[L]$ 이 최소화 하려면 0 이 되어야 한다.
   * 따라서 $\mathbb E[L] = \mathbb E[t \mid {\bf x}]$ 가 된다.
   * 이것은 단순히 이전에 유도한 결과와 같고, least square predictor 최적해는 바로 조건부 평균과 같다.

* 두번째 항은 바로 ${\bf x}$ 에서 평균된 분포 $t$ 의 분산(variance) 이다.
* 이것은 타겟데이터의 기본적인 변동에 대한 것이며, 노이즈로 해석이 가능하다.
* 두번째 항은 $y({\bf x})$ 와 관계가 없기 때문에 loss function 의 최소값에는 관계가 없다.

##### regression 문제를 푸는 3가지 접근 방법

   * 분류 문제 처럼, regression 에서도 적절한 확률분포를 찿아야만 한다. 그런 다음 최적의 결정을 하거나 또는 직접 decision 할 수 있도록 모델을 구성하여야 한다. regression 문제를 푸는데 역시 3가지 접근 방법을 제시한다.

   * $({\bf a})$ 결합확률 $p({\bf x}, t)$ 를 결정하는 추론 문제를 해결 하고, 다음 조건 확률 $p(t \mid {\bf x})$ 를 구하여, 다음 1.89 주어진 조건 평균을 구하기 위하여 maginalized 한다.

   * $({\bf b})$ 먼저 조건 확률분포 $p(t \mid {\bf x})$ 를 결정하기 위한 추론 문제를 해결 하고, 다음 1.89 를 구하기 위해서 marginalize 한다

   * $({\bf c})$ 학습데이타로 부터 직접 regression function $y({\bf x})$ 를 구한다.

   * 3 가지 접근 방법에 대한 상대적인 장점은 앞서 분류 문제에서 와 같다.

   * squired loss 는 regression 의 loss 함수를 구하는 유일한 방법은 아니다.
   * 사실, squired loss 가 나쁜 결과를 가져오는 경우에는 우린 더 현명한 접근방법을 발전 시켜야 한다.

##### multimodal 분포

   * 조건 확율 분포에서 발생하는 중요한 예제는 inverse problem 문제를 해결하는데서 나타나듯이 조건 확률 분포 $p(t \mid {\bf x})$ 가 멀티모달이라는 것이다 (분포에서 peak 가 여러가 나타나는 경우).
   * Minkowski loss 라고 불리는 squred loss 의 일반화 시킨 경우를 보면

$$\mathbb E[L_q] = \iint \left|y({\bf x})-t\right|^q p({\bf x}, t) d{\bf x}dt \qquad{(1.91)}$$

   * $q = 2$ 일때, expected loss 가 된다.
   * 함수 $\left |y - t \right |^q$ 가 q 에 관한 $y - t$ 에 대한 그래프가 그림 1.29 에 있다.

   * 그림 1.29 $L_q = \left |y - t\right |^q$ 의 q 에 따른 그래프

   * $\mathbb E[L_q]$ 를 최소는 $q = 2$ 에서 contitional mean 이고, $q = 1$ 와 $q \to 0$ 에서 conditional median 이다.

max_stars_repo_path: prml/Chap 1-5. Decision Theory.ipynb

max_issues_repo_path: prml/Chap 1-5. Decision Theory.ipynb

max_forks_repo_path: prml/Chap 1-5. Decision Theory.ipynb


text:
```python
%matplotlib inline
```

Graph attention network
==================================

**Authors:** `Hao Zhang <https://github.com/sufeidechabei/>`_, `Mufei Li
<https://github.com/mufeili>`_, `Minjie Wang
<https://jermainewang.github.io/>`_  `Zheng Zhang
<https://shanghai.nyu.edu/academics/faculty/directory/zheng-zhang>`_

In this tutorial, you learn about a graph attention network (GAT) and how it can be
implemented in PyTorch. You can also learn to visualize and understand what the attention
mechanism has learned.

The research described in the paper `Graph Convolutional Network (GCN) <https://arxiv.org/abs/1609.02907>`_,
indicates that combining local graph structure and node-level features yields
good performance on node classification tasks. However, the way GCN aggregates
is structure-dependent, which can hurt its generalizability.

One workaround is to simply average over all neighbor node features as described in
the research paper `GraphSAGE
<https://www-cs-faculty.stanford.edu/people/jure/pubs/graphsage-nips17.pdf>`_.
However, `Graph Attention Network <https://arxiv.org/abs/1710.10903>`_ proposes a
different type of aggregation. GAN uses weighting neighbor features with feature dependent and
structure-free normalization, in the style of attention.

Introducing attention to GCN
----------------------------

The key difference between GAT and GCN is how the information from the one-hop neighborhood is aggregated.

For GCN, a graph convolution operation produces the normalized sum of the node features of neighbors.

\begin{align}h_i^{(l+1)}=\sigma\left(\sum_{j\in \mathcal{N}(i)} {\frac{1}{c_{ij}} W^{(l)}h^{(l)}_j}\right)\end{align}

where $\mathcal{N}(i)$ is the set of its one-hop neighbors (to include
$v_i$ in the set, simply add a self-loop to each node),
$c_{ij}=\sqrt{|\mathcal{N}(i)|}\sqrt{|\mathcal{N}(j)|}$ is a
normalization constant based on graph structure, $\sigma$ is an
activation function (GCN uses ReLU), and $W^{(l)}$ is a shared
weight matrix for node-wise feature transformation. Another model proposed in
`GraphSAGE
<https://www-cs-faculty.stanford.edu/people/jure/pubs/graphsage-nips17.pdf>`_
employs the same update rule except that they set
$c_{ij}=|\mathcal{N}(i)|$.

GAT introduces the attention mechanism as a substitute for the statically
normalized convolution operation. Below are the equations to compute the node
embedding $h_i^{(l+1)}$ of layer $l+1$ from the embeddings of
layer $l$.

  :width: 450px
  :align: center

\begin{align}\begin{align}
  z_i^{(l)}&=W^{(l)}h_i^{(l)},&(1) \\
  e_{ij}^{(l)}&=\text{LeakyReLU}(\vec a^{(l)^T}(z_i^{(l)}||z_j^{(l)})),&(2)\\
  \alpha_{ij}^{(l)}&=\frac{\exp(e_{ij}^{(l)})}{\sum_{k\in \mathcal{N}(i)}^{}\exp(e_{ik}^{(l)})},&(3)\\
  h_i^{(l+1)}&=\sigma\left(\sum_{j\in \mathcal{N}(i)} {\alpha^{(l)}_{ij} z^{(l)}_j }\right),&(4)
  \end{align}\end{align}

Explanations:

* Equation (1) is a linear transformation of the lower layer embedding $h_i^{(l)}$
  and $W^{(l)}$ is its learnable weight matrix.
* Equation (2) computes a pair-wise *un-normalized* attention score between two neighbors.
  Here, it first concatenates the $z$ embeddings of the two nodes, where $||$
  denotes concatenation, then takes a dot product of it and a learnable weight vector
  $\vec a^{(l)}$, and applies a LeakyReLU in the end. This form of attention is
  usually called *additive attention*, contrast with the dot-product attention in the
  Transformer model.
* Equation (3) applies a softmax to normalize the attention scores on each node's
  incoming edges.
* Equation (4) is similar to GCN. The embeddings from neighbors are aggregated together,
  scaled by the attention scores.

There are other details from the paper, such as dropout and skip connections.
For the purpose of simplicity, those details are left out of this tutorial. To see more details,
download the `full example <https://github.com/dmlc/dgl/blob/master/examples/pytorch/gat/gat.py>`_.
In its essence, GAT is just a different aggregation function with attention
over features of neighbors, instead of a simple mean aggregation.

GAT in DGL
----------

To begin, you can get an overall impression about how a ``GATLayer`` module is
implemented in DGL. In this section, the four equations above are broken down
one at a time.

```python
import torch
import torch.nn as nn
import torch.nn.functional as F

class GATLayer(nn.Module):
    def __init__(self, g, in_dim, out_dim):
        super(GATLayer, self).__init__()
        self.g = g
        # equation (1)
        self.fc = nn.Linear(in_dim, out_dim, bias=False)
        # equation (2)
        self.attn_fc = nn.Linear(2 * out_dim, 1, bias=False)

    def edge_attention(self, edges):
        # edge UDF for equation (2)
        z2 = torch.cat([edges.src['z'], edges.dst['z']], dim=1)
        a = self.attn_fc(z2)
        return {'e': F.leaky_relu(a)}

    def message_func(self, edges):
        # message UDF for equation (3) & (4)
        return {'z': edges.src['z'], 'e': edges.data['e']}

    def reduce_func(self, nodes):
        # reduce UDF for equation (3) & (4)
        # equation (3)
        alpha = F.softmax(nodes.mailbox['e'], dim=1)
        # equation (4)
        h = torch.sum(alpha * nodes.mailbox['z'], dim=1)
        return {'h': h}

    def forward(self, h):
        # equation (1)
        z = self.fc(h)
        self.g.ndata['z'] = z
        # equation (2)
        self.g.apply_edges(self.edge_attention)
        # equation (3) & (4)
        self.g.update_all(self.message_func, self.reduce_func)
        return self.g.ndata.pop('h')
```

Equation (1)
^^^^^^^^^^^^

\begin{align}z_i^{(l)}=W^{(l)}h_i^{(l)},(1)\end{align}

The first one shows linear transformation. It's common and can be
easily implemented in Pytorch using ``torch.nn.Linear``.

Equation (2)
^^^^^^^^^^^^

\begin{align}e_{ij}^{(l)}=\text{LeakyReLU}(\vec a^{(l)^T}(z_i^{(l)}|z_j^{(l)})),(2)\end{align}

The un-normalized attention score $e_{ij}$ is calculated using the
embeddings of adjacent nodes $i$ and $j$. This suggests that the
attention scores can be viewed as edge data, which can be calculated by the
``apply_edges`` API. The argument to the ``apply_edges`` is an **Edge UDF**,
which is defined as below:

```python
def edge_attention(self, edges):
    # edge UDF for equation (2)
    z2 = torch.cat([edges.src['z'], edges.dst['z']], dim=1)
    a = self.attn_fc(z2)
    return {'e' : F.leaky_relu(a)}
```

Here, the dot product with the learnable weight vector $\vec{a^{(l)}}$
is implemented again using PyTorch's linear transformation ``attn_fc``. Note
that ``apply_edges`` will **batch** all the edge data in one tensor, so the
``cat``, ``attn_fc`` here are applied on all the edges in parallel.

Equation (3) & (4)
^^^^^^^^^^^^^^^^^^

\begin{align}\begin{align}
  \alpha_{ij}^{(l)}&=\frac{\exp(e_{ij}^{(l)})}{\sum_{k\in \mathcal{N}(i)}^{}\exp(e_{ik}^{(l)})},&(3)\\
  h_i^{(l+1)}&=\sigma\left(\sum_{j\in \mathcal{N}(i)} {\alpha^{(l)}_{ij} z^{(l)}_j }\right),&(4)
  \end{align}\end{align}

Similar to GCN, ``update_all`` API is used to trigger message passing on all
the nodes. The message function sends out two tensors: the transformed ``z``
embedding of the source node and the un-normalized attention score ``e`` on
each edge. The reduce function then performs two tasks:

* Normalize the attention scores using softmax (equation (3)).
* Aggregate neighbor embeddings weighted by the attention scores (equation(4)).

Both tasks first fetch data from the mailbox and then manipulate it on the
second dimension (``dim=1``), on which the messages are batched.

```python
def reduce_func(self, nodes):
    # reduce UDF for equation (3) & (4)
    # equation (3)
    alpha = F.softmax(nodes.mailbox['e'], dim=1)
    # equation (4)
    h = torch.sum(alpha * nodes.mailbox['z'], dim=1)
    return {'h' : h}
```

Multi-head attention
^^^^^^^^^^^^^^^^^^^^

Analogous to multiple channels in ConvNet, GAT introduces **multi-head
attention** to enrich the model capacity and to stabilize the learning
process. Each attention head has its own parameters and their outputs can be
merged in two ways:

\begin{align}\text{concatenation}: h^{(l+1)}_{i} =||_{k=1}^{K}\sigma\left(\sum_{j\in \mathcal{N}(i)}\alpha_{ij}^{k}W^{k}h^{(l)}_{j}\right)\end{align}

or

\begin{align}\text{average}: h_{i}^{(l+1)}=\sigma\left(\frac{1}{K}\sum_{k=1}^{K}\sum_{j\in\mathcal{N}(i)}\alpha_{ij}^{k}W^{k}h^{(l)}_{j}\right)\end{align}

where $K$ is the number of heads. You can use
concatenation for intermediary layers and average for the final layer.

Use the above defined single-head ``GATLayer`` as the building block
for the ``MultiHeadGATLayer`` below:

```python
class MultiHeadGATLayer(nn.Module):
    def __init__(self, g, in_dim, out_dim, num_heads, merge='cat'):
        super(MultiHeadGATLayer, self).__init__()
        self.heads = nn.ModuleList()
        for i in range(num_heads):
            self.heads.append(GATLayer(g, in_dim, out_dim))
        self.merge = merge

    def forward(self, h):
        head_outs = [attn_head(h) for attn_head in self.heads]
        if self.merge == 'cat':
            # concat on the output feature dimension (dim=1)
            return torch.cat(head_outs, dim=1)
        else:
            # merge using average
            return torch.mean(torch.stack(head_outs))
```

Put everything together
^^^^^^^^^^^^^^^^^^^^^^^

Now, you can define a two-layer GAT model.

```python
class GAT(nn.Module):
    def __init__(self, g, in_dim, hidden_dim, out_dim, num_heads):
        super(GAT, self).__init__()
        self.layer1 = MultiHeadGATLayer(g, in_dim, hidden_dim, num_heads)
        # Be aware that the input dimension is hidden_dim*num_heads since
        # multiple head outputs are concatenated together. Also, only
        # one attention head in the output layer.
        self.layer2 = MultiHeadGATLayer(g, hidden_dim * num_heads, out_dim, 1)

    def forward(self, h):
        h = self.layer1(h)
        h = F.elu(h)
        h = self.layer2(h)
        return h
```

We then load the Cora dataset using DGL's built-in data module.

```python
from dgl import DGLGraph
from dgl.data import citation_graph as citegrh
import networkx as nx

def load_cora_data():
    data = citegrh.load_cora()
    features = torch.FloatTensor(data.features)
    labels = torch.LongTensor(data.labels)
    mask = torch.ByteTensor(data.train_mask)
    g = data.graph
    # add self loop
    g.remove_edges_from(nx.selfloop_edges(g))
    g = DGLGraph(g)
    g.add_edges(g.nodes(), g.nodes())
    return g, features, labels, mask
```

The training loop is exactly the same as in the GCN tutorial.

```python
import time
import numpy as np

g, features, labels, mask = load_cora_data()

# create the model, 2 heads, each head has hidden size 8
net = GAT(g,
          in_dim=features.size()[1],
          hidden_dim=8,
          out_dim=7,
          num_heads=2)

# create optimizer
optimizer = torch.optim.Adam(net.parameters(), lr=1e-3)

# main loop
dur = []
for epoch in range(30):
    if epoch >= 3:
        t0 = time.time()

    logits = net(features)
    logp = F.log_softmax(logits, 1)
    loss = F.nll_loss(logp[mask], labels[mask])

    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

    if epoch >= 3:
        dur.append(time.time() - t0)

    print("Epoch {:05d} | Loss {:.4f} | Time(s) {:.4f}".format(
        epoch, loss.item(), np.mean(dur)))
```

Visualizing and understanding attention learned
----------------------------------------------

Cora
^^^^

The following table summarizes the model performance on Cora that is reported in
`the GAT paper <https://arxiv.org/pdf/1710.10903.pdf>`_ and obtained with DGL
implementations.

.. list-table::
   :header-rows: 1

   * - Model
     - Accuracy
   * - GCN (paper)
     - $81.4\pm 0.5%$
   * - GCN (dgl)
     - $82.05\pm 0.33%$
   * - GAT (paper)
     - $83.0\pm 0.7%$
   * - GAT (dgl)
     - $83.69\pm 0.529%$

*What kind of attention distribution has our model learned?*

Because the attention weight $a_{ij}$ is associated with edges, you can
visualize it by coloring edges. Below you can pick a subgraph of Cora and plot the
attention weights of the last ``GATLayer``. The nodes are colored according
to their labels, whereas the edges are colored according to the magnitude of
the attention weights, which can be referred with the colorbar on the right.

  :width: 600px
  :align: center

You can see that the model seems to learn different attention weights. To
understand the distribution more thoroughly, measure the `entropy
<https://en.wikipedia.org/wiki/Entropy_(information_theory>`_) of the
attention distribution. For any node $i$,
$\{\alpha_{ij}\}_{j\in\mathcal{N}(i)}$ forms a discrete probability
distribution over all its neighbors with the entropy given by

\begin{align}H({\alpha_{ij}}_{j\in\mathcal{N}(i)})=-\sum_{j\in\mathcal{N}(i)} \alpha_{ij}\log\alpha_{ij}\end{align}

A low entropy means a high degree of concentration, and vice
versa. An entropy of 0 means all attention is on one source node. The uniform
distribution has the highest entropy of $\log(\mathcal{N}(i))$.
Ideally, you want to see the model learns a distribution of lower entropy
(i.e, one or two neighbors are much more important than the others).

Note that since nodes can have different degrees, the maximum entropy will
also be different. Therefore, you plot the aggregated histogram of entropy
values of all nodes in the entire graph. Below are the attention histogram of
learned by each attention head.

|image2|

As a reference, here is the histogram if all the nodes have uniform attention weight distribution.

  :width: 250px
  :align: center

One can see that **the attention values learned is quite similar to uniform distribution**
(i.e, all neighbors are equally important). This partially
explains why the performance of GAT is close to that of GCN on Cora
(according to `author's reported result
<https://arxiv.org/pdf/1710.10903.pdf>`_, the accuracy difference averaged
over 100 runs is less than 2 percent). Attention does not matter
since it does not differentiate much.

*Does that mean the attention mechanism is not useful?* No! A different
dataset exhibits an entirely different pattern, as you can see next.

Protein-protein interaction (PPI) networks
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The PPI dataset used here consists of $24$ graphs corresponding to
different human tissues. Nodes can have up to $121$ kinds of labels, so
the label of node is represented as a binary tensor of size $121$. The
task is to predict node label.

Use $20$ graphs for training, $2$ for validation and $2$
for test. The average number of nodes per graph is $2372$. Each node
has $50$ features that are composed of positional gene sets, motif gene
sets, and immunological signatures. Critically, test graphs remain completely
unobserved during training, a setting called "inductive learning".

Compare the performance of GAT and GCN for $10$ random runs on this
task and use hyperparameter search on the validation set to find the best
model.

.. list-table::
   :header-rows: 1

   * - Model
     - F1 Score(micro)
   * - GAT
     - $0.975 \pm 0.006$
   * - GCN
     - $0.509 \pm 0.025$
   * - Paper
     - $0.973 \pm 0.002$

The table above is the result of this experiment, where you use micro `F1
score <https://en.wikipedia.org/wiki/F1_score>`_ to evaluate the model
performance.

Note

Below is the calculation process of F1 score:

  .. math::

     precision=\frac{\sum_{t=1}^{n}TP_{t}}{\sum_{t=1}^{n}(TP_{t} +FP_{t})}

     recall=\frac{\sum_{t=1}^{n}TP_{t}}{\sum_{t=1}^{n}(TP_{t} +FN_{t})}

     F1_{micro}=2\frac{precision*recall}{precision+recall}

  * $TP_{t}$ represents for number of nodes that both have and are predicted to have label $t$
  * $FP_{t}$ represents for number of nodes that do not have but are predicted to have label $t$
  * $FN_{t}$ represents for number of output classes labeled as $t$ but predicted as others.
  * $n$ is the number of labels, i.e. $121$ in our case.

During training, use ``BCEWithLogitsLoss`` as the loss function. The
learning curves of GAT and GCN are presented below; what is evident is the
dramatic performance adavantage of GAT over GCN.

  :width: 300px
  :align: center

As before, you can have a statistical understanding of the attentions learned
by showing the histogram plot for the node-wise attention entropy. Below are
the attention histograms learned by different attention layers.

*Attention learned in layer 1:*

|image5|

*Attention learned in layer 2:*

|image6|

*Attention learned in final layer:*

|image7|

Again, comparing with uniform distribution:

  :width: 250px
  :align: center

Clearly, **GAT does learn sharp attention weights**! There is a clear pattern
over the layers as well: **the attention gets sharper with a higher
layer**.

Unlike the Cora dataset where GAT's gain is minimal at best, for PPI there
is a significant performance gap between GAT and other GNN variants compared
in `the GAT paper <https://arxiv.org/pdf/1710.10903.pdf>`_ (at least 20 percent),
and the attention distributions between the two clearly differ. While this
deserves further research, one immediate conclusion is that GAT's advantage
lies perhaps more in its ability to handle a graph with more complex
neighborhood structure.

What's next?
------------

So far, you have seen how to use DGL to implement GAT. There are some
missing details such as dropout, skip connections, and hyper-parameter tuning,
which are practices that do not involve DGL-related concepts. For more information
check out the full example.

* See the optimized `full example <https://github.com/dmlc/dgl/blob/master/examples/pytorch/gat/gat.py>`_.
* The next tutorial describes how to speedup GAT models by parallelizing multiple attention heads and SPMV optimization.

.. |image2| image:: https://s3.us-east-2.amazonaws.com/dgl.ai/tutorial/gat/cora-attention-hist.png
.. |image5| image:: https://s3.us-east-2.amazonaws.com/dgl.ai/tutorial/gat/ppi-first-layer-hist.png
.. |image6| image:: https://s3.us-east-2.amazonaws.com/dgl.ai/tutorial/gat/ppi-second-layer-hist.png
.. |image7| image:: https://s3.us-east-2.amazonaws.com/dgl.ai/tutorial/gat/ppi-final-layer-hist.png

ext: ipynb


text:
# Lecture 5: Linear systems

## Previous lecture
- Matrix rank
- Skeleton decomposition
- Low-rank approximation
- Singular Value Decomposition (SVD)

## Today lecture
Today we talk about:
- Linear systems
- Inverse matrix
- Condition number
- Gaussian elimination

## Linear systems

Linear systems of equations are the basic tool in NLA.

They appear as:

- Linear regresion problems
- Discretization of partial differential/integral equations
- Linearizations of nonlinear regression problems
- Optimization (i.e., Gauss-Newton and Newton-Raphson methods, KKT).

## Linear equations and matrices
From school we know about linear equations.

A linear system of equations can be written in the form
\begin{equation}
\begin{split}
    &2 x + 3 y = 5\quad &\longrightarrow \quad &2x + 3 y + 0 z = 5\\
    &2 x + 3z = 5\quad &\longrightarrow\quad &2 x + 0 y + 3 z = 5\\
    &x + y = 2\quad &\longrightarrow\quad  & 1 x + 1 y + 0 z = 2\\
\end{split}
\end{equation}

## Matrix form
The "matrix form" is
$$
\begin{pmatrix}
2 & 3 & 0 \\
2 & 0 & 3 \\
1 & 1 & 0 \\
\end{pmatrix}\begin{pmatrix}
x \\
y \\
z
\end{pmatrix} =
\begin{pmatrix}
5 \\
5 \\
2
\end{pmatrix}
$$
or simply
$$
A u = f,
$$

where $A$ is a $3 \times 3$ matrix.

## Over/under determined linear systems

If the system $Au = f$ has more equations than unknowns it is called **overdetermined** system (generically, no solution)

If the system $Au = f$ has less equations than unknowns it is called **underdetermined** system (solution is non-unique, to make it unique additional assumptions have to be made).

## Existence of solutions

A solution to the linear system of equations with a **square** matrix $A$

$$A u = f$$

exists, iff
* $\det A \ne 0$

or

* Matrix $A$ has full rank.

## Scales of linear systems

In different applications, the typical size of the linear systems can be different.

- Small: $n \leq 10^4$ (full matrix can be stored in memory, **dense matrix**)
- Medium: $n = 10^4 - 10^6$ (typically, **sparse** or **structured**)
- Large: $n = 10^8 - 10^9$ (typicall **sparse** + parallel computations)

## Linear systems can be big
We take a continious problem, discretize on a mesh with $N$ elements and get a linear system with $N\times N$ system.
Example of a mesh around A319 aircraft
(taken from [GMSH website](http://geuz.org/gmsh/)).

The main difficulty is that these systems are big: millions or billions of unknowns!

## Linear systems can be structured
Storing $N^2$ elements of a matrix is prohibitive even for $N = 100000$.

How to work with such matrices?

Fortunately, those matrices are **structured** and require $\mathcal{O}(N)$ parameters to be stored.

The most widespread structure are **sparse matrices**: in the matrix there are $\mathcal{O}(N)$ non-zeros!

Example (one of the famous matrices around for $n = 5$):
$$
  \begin{pmatrix}
  2 & -1 & 0 & 0 & 0 \\
  -1 & 2 & -1 & 0 & 0 \\
  0 & -1 & 2 & -1 & 0 \\
  0 & 0 &-1& 2 & -1  \\
  0 & 0 & 0 & -1 & 2 \\
  \end{pmatrix}
$$
At least you can store such matrices, and multiply by vector fast; but how to solve linear systems?! (and that is typically the final goal).

## Main questions about linear systems

1. What is the accuracy we get from the solution (due to rounding errors)?
2. How we compute the solution? (LU-decomposition, Gaussian elimination)
3. What is the complexity of the solution of linear systems?

## How to solve linear systems?
**Important**: forget about determinants and the **Cramer rule** (it is  good for $2 \times 2$ matrices still)!

## How to solve linear systems?

The main tool is variable elimination.
\begin{equation}
\begin{split}
    &2 y + 3 x = 5 \quad&\longrightarrow \quad &y = 5/2 -  3/2 x \\
    &2 x + 3z = 5 \quad&\longrightarrow\quad &z = 5/3 - 2/3 x\\
    &x + y = 2 \quad&\longrightarrow\quad  & 5/2 + 5/3 - (3/2 + 2/3) x = 2,\\
\end{split}
\end{equation}
and that is how you find $x$ (and all previous ones).

This process is called **Gaussian elimination** and is one of the most widely used algorithms.

## Gaussian elimination
Gaussian elimination consists of two steps:
1. Forward step
2. Backward step

## Forward step
In the forward step, we eliminate $x_1$:

$$
   x_1 = f_1 - (a_{12} x_2 + \ldots + a_{1n} x_n)/a_{11},
$$

and then substitute this into the equations $2, \ldots, n$.

Then we eliminate $x_2$ and so on from the second equation.

The important thing is that the **pivots** (that we divide over) are not equal to $0$.

## Backward step
In the backward step, we solve equation for $x_n$, put it into the equation for $x_{n-1}$ and so on, until we
compute all $x_i, i=1,\ldots, n$.

## Gaussian elimination and LU decomposition
Gaussian elimination is the computation of one of the most important matrix decompositions: **LU-decomposition**.

**Definition**: LU-decomposition of the squared matrix $A$ is the representation

$$A =  LU,$$

where $L$ is **lower triangular** and $U$ is **upper triangular** matrix (i.e. elements strictly above the diagonal are zero, elements strictly below the diagonal are zero)

This factorization is **non-unique**, so it is typical to require that the matrix $L$ has ones on the diagonal.

**Main goal** of the LU decomposition is to solve linear systems, because

$$
    A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f,
$$

and this reduces to the solution of two linear systems **forward step**
$$
     L y = f,
$$
and **backward step**
$$
   U x = y.
$$

Does $LU$ decomposition always exist?

## Complexity of the Gaussian elimination/LU decomposition
Each elimination step requires $\mathcal{O}(n^2)$ operations.

Thus, the cost of the naive algorithm is $\mathcal{O}(n^3)$.

**Think a little bit**: Can Strassen help here?

## Block LU-decomposition

We can try to compute **block** version of LU-decomposition:

$$\begin{bmatrix} A_{11} & A_{12} \\
A_{21} & A_{22}
\end{bmatrix} = \begin{bmatrix} L_{11} & 0 \\
L_{21} & L_{22}
\end{bmatrix} \begin{bmatrix} U_{11} & U_{12} \\
0 & U_{22}
\end{bmatrix}
$$
There are two basic operations: compute LU-factorization of half-matrices + matrix-by-matrix product.

## Computing LU-decomposition
In many cases, computing LU-decomposition once is a good idea!

Once the decomposition is found (it is $\mathcal{O}(n^3)$), then solving linear systems with $L$ and $U$ costs only $\mathcal{O}(n^2)$ operations.

**Check:** Solving linear systems with triangular matrices is easy (why?). How we compute the $L$ and $U$ factors?

Not let us see how it works in numpy (for the Hilbert matrix, of course).

## Existence of the LU-decomposition
The LU-decomposition algorithm does not fail

if **we do not divide by zero** at every step of the Gaussian elimination.

When it is so, for which class of matrices?

It is true for **strictly regular matrices**.

## Strictly regular matrices and LU decomposition
A matrix $A$ is called strictly regular, if all of its **principal minors**

(i.e, submatrices in the first $k$ rows and $k$ columns) are non-singular.

In this case, there always exists a LU-decomposition.  The reverse is also true (check!).

**Proof**. If there is a LU-decomposition, then the matrix is strictly regular

This follows from the fact that to get a minor you multiply a corresponding submatrix of $L$ by a corresponding submatrix of $U$,

and they are non-singular (since non-singularity of triangular matrices is equivalent to the fact that their diagonal elements are not equal to zero).

The other way can be proven by **induction**. Suppose that we know that for all $(n-1) \times (n-1)$ matrices will non-singular minors LU-decomposition exists.

Then, consider the block form
$$
    A = \begin{bmatrix}
          a & c^{\top} \\
          b & D
    \end{bmatrix},
$$
where $D$ is $(n-1) \times (n-1)$.

Then we do "block elimination":

$$
     \begin{bmatrix}
     1 & 0 \\
     -z & I
     \end{bmatrix}
     \begin{bmatrix}
     a & c^{\top} \\
     b & D
     \end{bmatrix}=
     \begin{bmatrix}
     a & c^{\top} \\
     0 & A_1
     \end{bmatrix},
$$

where $z = \frac{b}{a}, \quad A_1 = D - \frac{1}{a} b c^{\top}$.
We can show that $A_1$ is also strictly regular, thus it has (by induction) the LU-decomposition:
$$
   A_1 = L_1 U_1,
$$
And that gives the $LU$ decomposition of the original matrix.

**Corollary:** If $L$ is unit triangular (ones on the diagonal), then $LU$-decomposition is unique.

*Proof:* Indeed, $L_1 U_1 = L_2 U_2$ means $L_2^{-1} L_1 = U_2 U_1^{-1}$. $L_2^{-1} L_1 $ is lower triangular with ones on the diagonal. $U_2 U_1^{-1}$ is upper triangular. Thus, $L_2^{-1} L_1 = U_2 U_1^{-1} = I$ and $L_1 = L_2$, $U_1 = U_2$.

## LU for positive definite Hermitian matrices (Cholesky factorization)

Strictly regular matrices have LU-decomposition.

An important **subclass** of strictly regular matrices is the class of **Hermitian positive definite matrices**

**Definition:** A matrix is called  positive definite  if for any $x: \Vert x \Vert \ne 0$ we have

$$(Ax, x) > 0.$$

**Claim:** A Hermitian positive definite matrix $A$ is strictly regular and has **Cholesky factorization** of the form

$$A = RR^*,$$

where $R$ is a lower triangular matrix.

Let us try to prove this fact (on the whiteboard).

It is sometimes referred to as "square root" of the matrix.

## When LU fails
What happens, if the matrix is not strictly regular (or the **pivots** in the Gaussian elimination are really small?).

There is classical $2 \times 2$ example of a matrix with a bad LU decomposition.

The matrix we look at is
$$
    A = \begin{pmatrix}
    \varepsilon & 1 \\
    1 & 1
    \end{pmatrix}
$$

If $\varepsilon$ is sufficiently small, we **might** fail. In contrast the Cholesky factorization is **always stable**.

Let us do some demo here.

```python
import numpy as np
eps = 1e-18#1.12e-16
a = [[eps, 1],[1.0,  1]]
a = np.array(a)
a0 = a.copy()
n = a.shape[0]
L = np.zeros((n, n))
U = np.zeros((n, n))
for k in range(n): #Eliminate one row
    L[k, k] = 1
    for i in range(k+1, n):
        L[i, k] = a[i, k] / a[k, k]
        for j in range(k+1, n):
            a[i, j] = a[i, j] - L[i, k] * a[k, j]
    for j in range(k, n):
        U[k, j] = a[k, j]

print 'L * U - A:', np.dot(L, U) - a0
L
```

    L * U - A: [[  0.00000000e+00   0.00000000e+00]
     [ -1.11022302e-16  -1.00000000e+00]]

    array([[  1.00000000e+00,   0.00000000e+00],
           [  1.00000000e+18,   1.00000000e+00]])

## The concept of pivoting
We can do pivoting, i.e. permute rows and columns to maximize $A_{kk}$ that we divide over.

The simplest but effective strategy is the **row pivoting**: at each step, select the index that is maximal in modulus, and put it onto the diagonal.

It gives us the decomposition

$$A = P L U,$$
where $P$ is a **permutation matrix**.

What makes **row pivoting** good?

It is made good by the fact that

$$
   | L_{ij}|<1,
$$
but the elements of $U$ can grow, up to $2^n$! (in practice, this is very rarely encountered).

Can you come up with a matrix where the elements of $U$ grow as much as possible?

# Stability of linear systems

There is a fundamental problem of solving linear systems which is independent on the algorithm used.

It occures when elements of a matrix are represented as floating point numbers or there is some measurement noise.

Let us illustrate this issue on the following example.

```python
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
n = 15
a = [[1.0/(i + j + 1) for i in xrange(n)] for j in xrange(n)]
a = np.array(a)
rhs = np.random.randn(n) #Right-hand side
x = np.linalg.solve(a, rhs) #This function computes LU-factorization and solves linear system

#And check if everything is fine
er = np.linalg.norm(a.dot(x) - rhs) / np.linalg.norm(rhs)
print er
plt.xkcd()
plt.plot(x)
```

As you see, the error grows with larger $n$, and we have to find out why.
**Important point** is that it is not a problem of the algorithm: it is a problem of representing
the matrix in the memory. The error occurs in the moment when the matrix elements are evaluated approximately.

## Linear systems and inverse matrix
What was the problem in the previous example?

Why the error grows so quickly?

And here is one of the main concepts of numerical linear algebra: the concept of **condition number** of a matrix.

But before that we have to define the **inverse**.

## Inverse: definition
The inverse of a matrix $A$ is defined as a matrix $X$ denoted by $A^{-1}$ such that
$$
   AX = XA = I,
$$
where $I$ is the identity matrix (i.e., $I_{ij} = 0$ if $i \ne j$ and $1$ otherwise).
The computation of the inverse is linked to the solution of linear systems.  Indeed, $i$-th column of the product gives
$$
A x_i = e_i,
$$
where $e_i$ is the $i$-th column of the identity matrix. Thus, we can apply Gaussian elimination to solve this system. Moreover, if there are no divisions by zero in this process (and the pivots do not depend on the right-hand side), then it is possible to solve the system.

## Inverse matrix and linear systems
If we have computed $A^{-1}$, the solution of linear system

$$Ax = f$$

is just $x = A^{-1} f$.

Indeed,

$$
    A(A^{-1} f) = (AA^{-1})f = I f = f.
$$

## Neumann series
To study, why there can be such big errors in a solution (see the example above on the Hilbert matrix)  we need an important auxilary result:

**Neumann series**:

If for a matrix $\Vert F \Vert < 1$ then the matrix $(I - F)$ is invertible and

$$(I - F)^{-1} = I + F + F^2 + F^3 + \ldots = \sum_{k=0}^{\infty} F^k.$$

Note that it is a matrix version of the geometric progression.

**Question**: What norm is considered here? What is the "best possible" norm here?

## Proof
The proof is constructive. First of all, the show prove that the series $\sum_{k=0}^{\infty} F^k$ converges.

Like in the scalar case, we have
$$
   (I - F) \sum_{k=0}^N F^k = (I - F^{N+1}) \rightarrow I.
$$

Indeed,
$$
\| (I - F^{N+1}) - I\| = \|F^{N+1}\| \leqslant \|F\|^{N+1} \to 0, \quad N\to +\infty.
$$

We can also estimate the **norm of the inverse**:
$$
  \Vert \sum_{k=0}^N F^k \Vert \leq \sum_{k=0}^N \Vert F \Vert^k \Vert I \Vert \leq \frac{\Vert I \Vert}{1 - \Vert F \Vert}
$$

## Small perturbation of the inverse
Using this result, we can estimate, how the perturbation of the matrix influences the inverse matrix. We assume
that the perturbation $E$ is small in the sense that $\Vert A^{-1} E \Vert < 1$. Then
$$(A + E)^{-1} = \sum_{k=0}^{\infty} (-A^{-1} E)^k A^{-1}$$
and moreover,
$$
  \frac{\Vert (A + E)^{-1} - A^{-1} \Vert}{\Vert A^{-1} \Vert} \leq \frac{\Vert A^{-1} \Vert \Vert E \Vert \Vert I \Vert}{1 - \Vert A^{-1} E \Vert}.
$$
As you see, the norm of the inverse enters the estimate.

## Condition number of a linear system
Now consider the **perturbed** linear system:
$$
   (A + \Delta A) \widehat{x} = f + \Delta f.
$$

## Estimates

$$
\begin{split}
\widehat{x} - x &= (A + \Delta A)^{-1} (f + \Delta f) - A^{-1} f =\\
&= \left((A + \Delta A)^{-1} - A^{-1}\right)f + (A + \Delta A)^{-1} \Delta f = \\
&= \Big[\sum_{k=0}^{\infty} (-A^{-1} \Delta A)^k\Big] A^{-1} f + \Big[\sum_{k=0}^{\infty} (A^{-1} \Delta A)^k \Big] A^{-1} \Delta f,
\end{split}
$$
therefore
$$
\begin{split}
\frac{\Vert \widehat{x} - x \Vert}{\Vert x \Vert} \leq
&\frac{\Vert A \Vert \Vert A^{-1} \Vert}{1 - \|A^{-1}\Delta A\|} \Big(\frac{\Vert\Delta A\Vert}{\Vert A \Vert} + \frac{\Vert \Delta f \Vert}{ \Vert f \Vert}\Big) \leq \\
\leq
&\frac{\Vert A \Vert \Vert A^{-1} \Vert}{1 - \|A\|\|A^{-1}\|\frac{\|\Delta A\|}{\|A\|}} \Big(\frac{\Vert\Delta A\Vert}{\Vert A \Vert} + \frac{\Vert \Delta f \Vert}{ \Vert f \Vert}\Big) \equiv \\
\equiv &\frac{\mathrm{cond}(A)}{1 - \mathrm{cond}(A)\frac{\|\Delta A\|}{\|A\|}} \Big(\frac{\Vert\Delta A\Vert}{\Vert A \Vert} + \frac{\Vert \Delta f \Vert}{ \Vert f \Vert}\Big)
\end{split}
$$

The crucial role is played by the **condition number** $\mathrm{cond}(A) = \Vert A \Vert \Vert A^{-1} \Vert$.

## Condition number
The larger the condition number, the less number of digits we can recover. Note, that the condition number is different for different norms.

Note, that if $\Delta A = 0$, then
$$
\frac{\Vert \widehat{x} - x \Vert}{\Vert x \Vert} \leq \mathrm{cond}(A) \frac{\|\Delta f\|}{\|f\|}
$$

The spectral norm of the matrix is equal to **largest singular value**, and the singular values of the inverse matrix are equal to the inverses of the singular values.  Thus,  the condition number is equal to the ratio of the largest singular value and the smallest singular value.
$$
    \mathrm{cond}_2 (A) = \|A\|_2 \|A^{-1}\|_2 = \frac{\sigma_{\max}}{\sigma_{\min}}
$$

## Hilbert matrix (again)
We can also try to test how tight is the estimate, both with ones in the right-hand side, and with a random vector in the right-hand side. The results are strickingly different

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

n = 1000
a = [[1.0/(i + j + 1) for i in xrange(n)] for j in xrange(n)]
a = np.array(a)
rhs = np.ones(n) #Right-hand side
f = np.linalg.solve(a, rhs)

#And check if everything is fine
er = np.linalg.norm(a.dot(f) - rhs) / np.linalg.norm(rhs)
cn = np.linalg.cond(a, 2)
print 'Error:', er, 'Condition number:', cn
```

    Error: 8.96141017158e-08 Condition number: 1.5912638327e+20

And with random right-hand side...

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

n = 100
a = [[1.0/(i + j + 1) for i in xrange(n)] for j in xrange(n)]
a = np.array(a)
rhs = np.random.randn(n) #Right-hand side
f = np.linalg.solve(a, rhs)

#And check if everything is fine
er = np.linalg.norm(a.dot(f) - rhs) / np.linalg.norm(rhs)
cn = np.linalg.cond(a)
print 'Error:', er, 'Condition number:', cn

u, s, v = np.linalg.svd(a)
rhs = np.random.randn(n)
plt.plot(u.T.dot(rhs))
```

** Can you think about an explanation?**

## Overdetermined linear systems
Important class of problems are **overdetermined linear systems**, when the number of equations is greater, than the number of unknowns. The simplest example that you all know, is **linear fitting**, fitting a set of 2D points by a line.

Then, a typical way is to minimize the residual (**least squares**)

$$\Vert A x - b \Vert_2 \rightarrow \min$$

## Overdetermined system and Gram matrix

The optimality condition is $0\equiv \delta \left(\|Ax-b\|_2^2\right)$, where $\delta$ denotes variation. Therefore,

$$
\begin{split}
0\equiv \delta \left(\|Ax-b\|_2^2\right) =\delta \left(Ax-b, Ax - b\right) = (A\delta x, Ax-b)+ (Ax-b, A\delta x) =
2(A x - b, A \delta x) = 2(A^*(A x - b), \delta x) = 0,
\end{split}
$$
for any $\delta x$. Here we used, that $\delta(Ax-b)=A\delta x$ ($A$ is linear operator and $b$ is constant) and $(x,Ay) \equiv (A^*x, y)$. Thus,
$$
\quad A^* A x = A^* b
$$
The matrix $A^* A$ is called **Gram matrix** and the system is called **normal equations**.

This is not a good way to do it, since the condition number of $A^* A$ is a square of condition number of $A$ (check why).

## Pseudoinverse

Matrix $A^* A$ can be singular in general case.
Therefore, we need to introduce the concept of pseudoinverse matrix  $A^{\dagger}$ such that

solution to the linear least squares problem can formally be written as
$$x = A^{\dagger} b.$$

The matrix $$A^{\dagger} = \lim_{\alpha \rightarrow 0}(\alpha I + A^* A)^{-1} A^*$$ is called **Moore-Penrose** pseudoinverse of the matrix $A$.

* If matrix $A$ has full column rank, then $A^* A$ is non-singular and we get $A^{\dagger} = \lim_{\alpha \rightarrow 0}(\alpha I + A^* A)^{-1} A^* = (A^* A)^{-1} A^*$.

* If matrix $A$ is squared and non-singular we get $A^{\dagger} = \lim_{\alpha \rightarrow 0}(\alpha I + A^* A)^{-1} A^* = (A^* A)^{-1} A^* = A^{-1} A^{-*} A^* = A^{-1}$ - standard inverse of $A$

* If $A$ has linearly dependent columns, then $A^\dagger b$ gives solution that has minimal Euclidean norm

## Compute pseudoinverse via SVD
Let $A = U \Sigma V^*$ be the SVD of $A$. Then,

$$A^{\dagger} = V \Sigma^{\dagger} U^*,$$

where $\Sigma^{\dagger}$ consists of inverses of non-zero singular values of $A$. Indeed,

$$A^{\dagger} = \lim_{\alpha \rightarrow 0}(\alpha I + A^* A)^{-1} A^* = \lim_{\alpha \rightarrow 0}( \alpha VV^* + V \Sigma^2 V^*)^{-1} V \Sigma U^* = \lim_{\alpha \rightarrow 0}( V(\alpha I + \Sigma^2) V^*)^{-1} V \Sigma U^* = V \lim_{\alpha \rightarrow 0}(\alpha I + \Sigma^2)^{-1} \Sigma U^* = V \Sigma^{\dagger} U^*,$$

* One can check that $\Sigma^{\dagger}$ contains just the inversion of nonzero singular values.

* If singular values are small one can skip inverting them. This will result in a solution which is less sensitive to the noise in the right-hand-side.
* The condition number for the Euclidean norm is still just the ratio of largest and smallest non-zero singular values.

## A canonical way to solve linear least squares

Is to use the $QR$ decomposition.

Any matrix can be factored into a product

$$A = Q R, $$
where $Q$ is unitary, and $R$ is upper triangular (details in the next lectures).

Then, since $$\|Ax - b\|_2 = \|QRx - b\|_2 =\|Rx - Q^*b\|_2$$ finding optimal $x$ is equivalent to solving $Rx = Q^* b$.

## Padding into a bigger system

Instead of solving $A^* A x = A^* b$,

we introduce a new variable $r = Ax - b$ and then have

$$A^* r = 0, \quad r = Ax - b,$$

or in the block form

$$
\begin{bmatrix}
  0 & A^* \\
  A & -I
\end{bmatrix}
\begin{bmatrix}
x \\
r
\end{bmatrix}
=
\begin{bmatrix}
0 \\
b
\end{bmatrix},
$$

the total size of the system is $(n + m)$ square, and the condition number is the same as for $A$ (by the way, how we define the condition number of a rectangular matrix?)

## Example of LS
Consider a two-dimensional example. Suppose we have a linear model
$$y = ax + b$$
and noisy data $(x_1, y_1), \dots (x_n, y_n)$. Then the linear system on coefficients will look as follows
$$
\begin{split}
a x_1 &+ b &= y_1 \\
&\vdots \\
a x_n &+ b &= y_n \\
\end{split}
$$
or in a matrix form
$$
\begin{pmatrix}
x_1 & 1 \\
\vdots & \vdots \\
x_n & 1 \\
\end{pmatrix}
\begin{pmatrix}
a \\
b
\end{pmatrix} =
\begin{pmatrix}
y_1 \\
\vdots  \\
y_n \\
\end{pmatrix},
$$
which represents overdetermined system.

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

a_exact = 1.
b_exact = 2.

n = 10
xi = np.arange(n)
yi = a_exact * xi + b_exact + 2*np.random.random(n)

A = np.array([xi, np.ones(n)])
coef = np.linalg.pinv(A).T.dot(yi) # coef is [a, b]

plt.xkcd()
plt.plot(xi, yi, 'o', label='$(x_i, y_i)$')
plt.plot(xi, coef[0]*xi + coef[1], label='Least Squares')
plt.legend(loc='best')
```

## Lacking for structure
A typical 3D-problem requires a $100 \times 100 \times 100$ discretization

This gives a linear system with $10^6$ unknowns, right hand side takes $8$ megabytes of memory.

This matrix has $10^6 \times 10^6 = 10^{12}$ elements, takes $8$ terabytes of memory.

Fortunately, the matrices in real-life are not **dense**, but have certain **structure**:

- Sparse matrices
- Low-rank matrices
- Toeplitz matrices (shift-invariant property)
- Sparse in certain bases

## Summary
- Linear systems can be solved by Gaussian elimination, complexity in $\mathcal{O}(n^3)$.
- Linear systems can be solved by LU-decomposition, complexity is $\mathcal{O}(n^3)$ for the decomposition, $\mathcal{O}(n^2)$ for each solve
- Linear least squares can be solved by normal equation (bad)
- Linear least squares can be solved by QR-decomposition (good) or by augmentation (not bad)
- Without structure, we can solve up to $10^4$ linear systems on a laptop (memory restrictions)

## Next lecture
- Eigenvectors, eigenvalues, Schur theorem

##### Questions?

```python
from IPython.core.display import HTML
def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>


text:
# For Three-pendulum, the Hamiltonian is

Three-pendulum system has three masses, $m_1, m_2, m_3$, three strings $L_1, L_2, L_3$ and three angles $\phi_1, \phi_2, \phi_3$

Through calculation for each mass's poential energy and kinetic energy, the $\mathcal{H}$ is simply $T + V$.

Hamilton's equations are

$\begin{align}
  \dot\phi_1 &=  \frac{p_\phi(1)}{m_1L_1} \\
  \dot\phi_2 &=  \frac{p_\phi(2)}{m_2L_2} \\
  \dot\phi_3 &=  \frac{p_\phi(3)}{m_3L_3} \\
  \dot p_\phi(1) &= -\frac{\partial\mathcal{H}}{\partial \phi_1} = -(m_1+m_2+m_3)gL_1sin(\phi_1)+(m_2+m_3)L_1L_2\dot\phi_1\dot\phi_2-m_3L_1L_2L_3\dot\phi_1\dot\phi_2\dot\phi_3sin(\phi_3-\phi_1+\phi_2)
  \\
   \dot p_\phi(2) &= -\frac{\partial\mathcal{H}}{\partial \phi_2} = -(m_2+m_3)gL_2sin(\phi_2)-(m_2+m_3)L_1L_2\dot\phi_1\dot\phi_2+m_3L_1L_2L_3\dot\phi_1\dot\phi_2\dot\phi_3sin(\phi_3-\phi_1+\phi_2)
  \\
   \dot p_\phi(3) &= -\frac{\partial\mathcal{H}}{\partial \phi_3} = -m_3gL_3sin(\phi_3)+m_3L_1L_2L_3\dot\phi_1\dot\phi_2\dot\phi_3sin(\phi_3-\phi_1+\phi_2)
  \;.
\end{align}$

```python
%matplotlib inline
```

```python
import numpy as np
from scipy.integrate import solve_ivp

import matplotlib.pyplot as plt

# Import explicitly the widgets we might use
import ipywidgets as widgets
from ipywidgets import HBox, VBox, Layout, Tab, Label, Checkbox
from ipywidgets import FloatSlider, IntSlider, Play, Dropdown, HTMLMath
from IPython.display import display

# Animation
from matplotlib import animation, rc
from IPython.display import HTML
```

```python
# The dpi (dots-per-inch) setting will affect the resolution and how large
#  the plots appear on screen and printed.  So you may want/need to adjust
#  the figsize when creating the figure.
plt.rcParams['figure.dpi'] = 100.    # this is the default for notebook

# Change the common font size (smaller when higher dpi)
font_size = 10
plt.rcParams.update({'font.size': font_size})
```

## Pendulum class and utility functions

```python
class HamiltonianPendulum():
    """
    Pendulum class implements the parameters and Hamilton's equations for
     a simple pendulum (no driving or damping).

    Parameters
    ----------
        Length : float
        The length of the strings
        mass : float
        mass of pendulum

    Methods
    -------
    dy_dt(t, y)
        Returns the right side of the differential equation in vector y,
        given time t and the corresponding value of y.
    """
    def __init__(self, L1=1., L2=1., L3=1., m1=1., m2=1., m3=1., g=1.):
        self.L1 = L1
        self.L2 = L2
        self.L3 = L3
        self.m1 = m1
        self.m2 = m2
        self.m3 = m3
        self.g = g
    def dy_dt(self, t, y):
        """
        This function returns the right-hand sides of Hamilton's eqs
        [dq1/dt dq2/dt dq3/dt dp1/dt dp2/dt dp3/dt] --> [dphi1/dt dphi2/dt dphi3/dt dp_phi1/dt dp_phi2/dt dp_phi3/dt]

        Parameters
        ----------
        t : float
            time
        y : float
            y[0] = phi1(t) y[1] = phi2(t) y[2]= phi3(t) y[3]= p_phi1(t) y[4]=p_phi2(t) y[5]=p_phi3(t)

        Returns
        -------

        """
        return [y[3]/(self.m1*self.L1), y[4]/(self.m2*self.L2), y[5]/(self.m3*self.L3),\
                -( self.m1 + self.m2 + self.m3 ) * self.L1 * self.g * np.sin(y[0]) + (self.m2 + self.m3) * self.L1 * self.L2 * y[3]/self.m1 * y[4]/self.m2 * np.sin(y[0]-y[1]) - self.m3 * self.L1 * self.L2 * self.L3 * y[3]/self.m1 * y[4]/self.m2 * y[5]/self.m3 * np.sin(y[2] - y[0] + y[1]),\
                -( self.m2 + self.m3 ) * self.L2 * self.g * np.sin(y[1]) - (self.m2 + self.m3) * self.L1 * self.L2 * y[3]/self.m1 * y[4]/self.m2 * np.sin(y[0]-y[1]) + self.m3 * self.L1 * self.L2 * self.L3 * y[3]/self.m1 * y[4]/self.m2 * y[5]/self.m3 * np.sin(y[2] - y[0] + y[1]),\
                - self.m3 * self.L3 * self.g * np.sin(y[2]) + self.m3 * self.L1 * self.L2 * self.L3 * y[3]/self.m1 * y[4]/self.m2 * y[5]/self.m3 * np.sin(y[2] - y[0] + y[1])]

    def solve_ode(self, t_pts, q1_0, q2_0, q3_0, p1_0, p2_0, p3_0,
                  abserr=1.0e-10, relerr=1.0e-10):
        """
        Solve the ODE given initial conditions.
        Specify smaller abserr and relerr to get more precision.
        """
        y = [q1_0, q2_0, q3_0, p1_0, p2_0, p3_0]
        solution = solve_ivp(self.dy_dt, (t_pts[0], t_pts[-1]),
                             y, t_eval=t_pts,
                             atol=abserr, rtol=relerr)
        q1, q2, q3, p1, p2, p3 = solution.y

        return q1, q2, q3, p1, p2, p3
```

```python
def plot_y_vs_x(x, y, axis_labels=None, label=None, title=None,
                color=None, linestyle=None, semilogy=False, loglog=False,
                ax=None):
    """
    Generic plotting function: return a figure axis with a plot of y vs. x,
    with line color and style, title, axis labels, and line label
    """
    if ax is None:        # if the axis object doesn't exist, make one
        ax = plt.gca()

    if (semilogy):
        line, = ax.semilogy(x, y, label=label,
                            color=color, linestyle=linestyle)
    elif (loglog):
        line, = ax.loglog(x, y, label=label,
                          color=color, linestyle=linestyle)
    else:
        line, = ax.plot(x, y, label=label,
                    color=color, linestyle=linestyle)

    if label is not None:    # if a label if passed, show the legend
        ax.legend()
    if title is not None:    # set a title if one if passed
        ax.set_title(title)
    if axis_labels is not None:  # set x-axis and y-axis labels if passed
        ax.set_xlabel(axis_labels[0])
        ax.set_ylabel(axis_labels[1])

    return ax, line
```

```python
def start_stop_indices(t_pts, plot_start, plot_stop):
    start_index = (np.fabs(t_pts-plot_start)).argmin()  # index in t_pts array
    stop_index = (np.fabs(t_pts-plot_stop)).argmin()  # index in t_pts array
    return start_index, stop_index
```

```python

```

## Make simple pendulum plots

```python
# Common plotting time (generate the full time then use slices)
t_start = 0.
t_end = 50.
delta_t = 0.001

t_pts = np.arange(t_start, t_end+delta_t, delta_t)

L1 = 1.
L2 = 1.
L3 = 1.
g = 1.
m1 = 1.
m2 = 1.
m3 = 1.
# Instantiate a pendulum
p1 = HamiltonianPendulum(L1=L1, L2=L2, L3=L3, m1=m1, m2=m2, m3=m3, g=g)

```

```python
# same initial conditions
phi1_0 = 0.5
phi2_0 = 0.5
phi3_0 = 0.5
p_phi1_0 = 0.
p_phi2_0 = 0.
p_phi3_0 = 0.

phi1, phi2, phi3, p_phi1, p_phi2, p_phi3  = p1.solve_ode(t_pts, phi1_0, phi2_0, phi3_0, p_phi1_0, p_phi2_0, p_phi3_0)

# start the plot!
fig = plt.figure(figsize=(15,5))

# first plot: phi1 plot
ax_a = fig.add_subplot(2,3,1)

start, stop = start_stop_indices(t_pts, t_start, t_end)
plot_y_vs_x(t_pts[start : stop], phi1[start : stop],
            color='blue',
            label=None,
            title=r'$\phi1(t)$',
            ax=ax_a)

# Second plot: phi2 plot
ax_b = fig.add_subplot(2,3,2)

start, stop = start_stop_indices(t_pts, t_start, t_end)
plot_y_vs_x(t_pts[start : stop], phi2[start : stop],
            color='blue',
            label=None,
            title=r'$\phi2(t)$',
            ax=ax_b)

# third plot: phi3 plot
ax_c = fig.add_subplot(2,3,3)

start, stop = start_stop_indices(t_pts, t_start, t_end)
plot_y_vs_x(t_pts[start : stop], phi3[start : stop],
            color='blue',
            label=None,
            title=r'$\phi3(t)$',
            ax=ax_c)

# Fourth plot: p_phi1 plot
ax_d = fig.add_subplot(2,3,4)

start, stop = start_stop_indices(t_pts, t_start, t_end)
plot_y_vs_x(t_pts[start : stop], p_phi1[start : stop],
            color='blue',
            label=None,
            title=r'$p1(t)$',
            ax=ax_d)

# Fifth plot: p_phi2 plot
ax_e = fig.add_subplot(2,3,5)

start, stop = start_stop_indices(t_pts, t_start, t_end)
plot_y_vs_x(t_pts[start : stop], p_phi2[start : stop],
            color='blue',
            label=None,
            title=r'$p2(t)$',
            ax=ax_e)

# Sixth plot: p_phi3 plot
ax_f = fig.add_subplot(2,3,6)

start, stop = start_stop_indices(t_pts, t_start, t_end)
plot_y_vs_x(t_pts[start : stop], p_phi3[start : stop],
            color='blue',
            label=None,
            title=r'$p3(t)$',
            ax=ax_f)

fig.tight_layout()
fig.savefig('simple_pendulum_Hamilton.png', bbox_inches='tight')
```

Now trying the power spectrum, plotting only positive frequencies and cutting off the lower peaks:

```python
start, stop = start_stop_indices(t_pts, t_start, t_end)
signal = phi1[start:stop]
power_spectrum = np.abs(np.fft.fft(signal))**2
freqs = np.fft.fftfreq(signal.size, delta_t)
idx = np.argsort(freqs)

fig_ps = plt.figure(figsize=(5,5))
ax_ps = fig_ps.add_subplot(1,1,1)
ax_ps.semilogy(freqs[idx], power_spectrum[idx], color='blue')
ax_ps.set_xlim(0, 1.)
ax_ps.set_ylim(1.e3, 1.e12)
ax_ps.set_xlabel('frequency')
ax_ps.set_title('Power Spectrum')
ax_ps.axvline(1./(2.*np.pi))
```

```python
start, stop = start_stop_indices(t_pts, t_start, t_end)
signal = phi2[start:stop]
power_spectrum = np.abs(np.fft.fft(signal))**2
freqs = np.fft.fftfreq(signal.size, delta_t)
idx = np.argsort(freqs)

fig_ps = plt.figure(figsize=(5,5))
ax_ps = fig_ps.add_subplot(1,1,1)
ax_ps.semilogy(freqs[idx], power_spectrum[idx], color='blue')
ax_ps.set_xlim(0, 1.)
ax_ps.set_ylim(1.e3, 1.e12)
ax_ps.set_xlabel('frequency')
ax_ps.set_title('Power Spectrum')
ax_ps.axvline(1./(2.*np.pi))
```

```python
start, stop = start_stop_indices(t_pts, t_start, t_end)
signal = phi3[start:stop]
power_spectrum = np.abs(np.fft.fft(signal))**2
freqs = np.fft.fftfreq(signal.size, delta_t)
idx = np.argsort(freqs)

fig_ps = plt.figure(figsize=(5,5))
ax_ps = fig_ps.add_subplot(1,1,1)
ax_ps.semilogy(freqs[idx], power_spectrum[idx], color='blue')
ax_ps.set_xlim(0, 1.)
ax_ps.set_ylim(1.e3, 1.e12)
ax_ps.set_xlabel('frequency')
ax_ps.set_title('Power Spectrum')
ax_ps.axvline(1./(2.*np.pi))
```

```python
#Get the widget
def update_plot(phi1_0, phi2_0,phi3_0, p_phi1_0, p_phi2_0,p_phi3_0):
    phi1, phi2, phi3, p_phi1, p_phi2, p_phi3  = p1.solve_ode(t_pts, phi1_0, phi2_0, phi3_0, p_phi1_0, p_phi2_0, p_phi3_0)
    # start the plot!
    fig = plt.figure(figsize=(15,5))
    # first plot: phi1 plot
    ax_a = fig.add_subplot(1,6,1)

    start, stop = start_stop_indices(t_pts, t_start, t_end)
    plot_y_vs_x(t_pts[start : stop], phi1[start : stop],
                color='blue',
                label=None,
                title=r'$\phi1(t)$',
                ax=ax_a)

    # Second plot: phi2 plot
    ax_b = fig.add_subplot(1,6,2)

    start, stop = start_stop_indices(t_pts, t_start, t_end)
    plot_y_vs_x(t_pts[start : stop], phi2[start : stop],
                color='blue',
                label=None,
                title=r'$\phi2(t)$',
                ax=ax_b)

    # third plot: phi3 plot
    ax_c = fig.add_subplot(1,6,3)

    start, stop = start_stop_indices(t_pts, t_start, t_end)
    plot_y_vs_x(t_pts[start : stop], phi3[start : stop],
                color='blue',
                label=None,
                title=r'$\phi3(t)$',
                ax=ax_c)

    # Fourth plot: p_phi1 plot
    ax_d = fig.add_subplot(1,6,4)

    start, stop = start_stop_indices(t_pts, t_start, t_end)
    plot_y_vs_x(t_pts[start : stop], p_phi1[start : stop],
                color='blue',
                label=None,
                title=r'$p1(t)$',
                ax=ax_d)

    # Fifth plot: p_phi2 plot
    ax_e = fig.add_subplot(1,6,5)

    start, stop = start_stop_indices(t_pts, t_start, t_end)
    plot_y_vs_x(t_pts[start : stop], p_phi2[start : stop],
                color='blue',
                label=None,
                title=r'$p2(t)$',
                ax=ax_e)

    # Sixth plot: p_phi3 plot
    ax_f = fig.add_subplot(1,6,6)

    start, stop = start_stop_indices(t_pts, t_start, t_end)
    plot_y_vs_x(t_pts[start : stop], p_phi3[start : stop],
                color='blue',
                label=None,
                title=r'$p3(t)$',
                ax=ax_f)

phi1_0_w = FloatSlider(value=1.0, min=0.1, max=10., step=0.1,
                      description=r'$\phi1_0_w$ :')
phi2_0_w = FloatSlider(value=1.0, min=0.1, max=10., step=0.1,
                      description=r'$\phi2_0_w$ :')
phi3_0_w = FloatSlider(value=1.0, min=0.1, max=10., step=0.1,
                      description=r'$\phi2_0_w$ :')
p_phi1_0_w = FloatSlider(value=1.0, min=0., max=10., step=0.1,
                      description=r'$\p_phi1_0_w$ :')
p_phi2_0_w = FloatSlider(value=1.0, min=0., max=10., step=0.1,
                      description=r'$\p_phi2_0_w$ :')
p_phi3_0_w = FloatSlider(value=1.0, min=0., max=10., step=0.1,
                      description=r'$\p_phi2_0_w$ :')
plot_out = widgets.interactive_output(update_plot,
                                      dict(
                                      phi1_0=phi1_0_w,
                                      phi2_0=phi2_0_w,
                                      phi3_0=phi3_0_w,
                                      p_phi1_0=p_phi1_0_w,
                                      p_phi2_0=p_phi2_0_w,
                                      p_phi3_0=p_phi3_0_w,

                                      )
                                     )

hbox_1 = HBox([phi1_0_w, phi2_0_w, phi3_0_w])
hbox_2 = HBox([p_phi1_0_w, p_phi2_0_w, p_phi3_0_w])
UI_box = VBox([hbox_1, hbox_2, plot_out])
display(UI_box)
```

    VBox(children=(HBox(children=(FloatSlider(value=1.0, description='$\\phi1_0_w$ :', max=10.0, min=0.1), FloatSl…

```python

```

```python

```

```python

```

```python

```


text:
.. _nb_custom:

## Custom Variable Type

In the following, we describe a custom variable problem. The variable is a string with a fixed length in our case.
However, we formulate the problem to be easily extended to have a variable length.
The objective function looks as follows:

\begin{align}
\begin{split}
\max f_1(x) & = & \, \# a \\[2mm]
\max f_2(x) & = & \, \# b
\end{split}
\end{align}

The first objective is the number of a's in the string and the second the number of b's.
For instance, for the variable "abdfgdgabb" the $f_1(x)=2$ and $f_2(x)=3$.

```python
import numpy as np
from pymoo.model.problem import Problem

class MyProblem(Problem):

    def __init__(self, n_characters=10):
        super().__init__(n_var=1, n_obj=2, n_constr=0, elementwise_evaluation=True)
        self.n_characters = n_characters
        self.ALPHABET = [c for c in string.ascii_lowercase]

    def _evaluate(self, x, out, *args, **kwargs):
        n_a, n_b = 0, 0
        for c in x[0]:
            if c == 'a':
                n_a += 1
            elif c == 'b':
                n_b += 1

        out["F"] = np.array([- n_a, - n_b], dtype=np.float)

```

The problem definition above define a problem with just one variable. This variable can be considered as a complex object which is in our case a string. The same principle can be used to use other data structure such as trees or lists with variable lengths. Because both objectives have to be maximized, we are minimizing their negative values.

To solve the problem the evolutionary operators sampling, crossover, mutation and duplication check needs to be implemented.
Each of the modules will be shown in the following.

### Sampling

Our sampling method just generates a random string which is equivalent of choosing a random letter from the alphabet (only lower case).
Because of the implementation of having only one variable, we return a matrix with the shape (n,1).

```python
from pymoo.model.sampling import Sampling

class MySampling(Sampling):

    def _do(self, problem, n_samples, **kwargs):
        X = np.full((n_samples, 1), None, dtype=np.object)

        for i in range(n_samples):
            X[i, 0] = "".join([np.random.choice(problem.ALPHABET) for _ in range(problem.n_characters)])

        return X

```

### Crossover

The crossover operator combines parents to create offsprings. In our framework the crossover operator retrieves the input already with predefined matings.
Our crossover randomly picks a character from the first or from the second parent.

```python
from pymoo.model.crossover import Crossover

class MyCrossover(Crossover):
    def __init__(self):

        # define the crossover: number of parents and number of offsprings
        super().__init__(2, 2)

    def _do(self, problem, X, **kwargs):

        # The input of has the following shape (n_parents, n_matings, n_var)
        _, n_matings, n_var = X.shape

        # The output owith the shape (n_offsprings, n_matings, n_var)
        # Because there the number of parents and offsprings are equal it keeps the shape of X
        Y = np.full_like(X, None, dtype=np.object)

        # for each mating provided
        for k in range(n_matings):

            # get the first and the second parent
            a, b = X[0, k, 0], X[1, k, 0]

            # prepare the offsprings
            off_a = ["_"] * problem.n_characters
            off_b = ["_"] * problem.n_characters

            for i in range(problem.n_characters):
                if np.random.random() < 0.5:
                    off_a[i] = a[i]
                    off_b[i] = b[i]
                else:
                    off_a[i] = b[i]
                    off_b[i] = a[i]

            # join the character list and set the output
            Y[0, k, 0], Y[1, k, 0] = "".join(off_a), "".join(off_b)

        return Y
```

### Mutation

The mutation needs to be implemented for our string object as well. We either change the order of the string (40%), randomly pick a new character with a given probability (40%) or leave the string unmodified (20%).

```python
from pymoo.model.mutation import Mutation

class MyMutation(Mutation):
    def __init__(self):
        super().__init__()

    def _do(self, problem, X, **kwargs):

        # for each individual
        for i in range(len(X)):

            r = np.random.random()

            # with a probabilty of 40% - change the order of characters
            if r < 0.4:
                perm = np.random.permutation(problem.n_characters)
                X[i, 0] = "".join(np.array([e for e in X[i, 0]])[perm])

            # also with a probabilty of 40% - change a character randomly
            elif r < 0.8:
                prob = 1 / problem.n_characters
                mut = [c if np.random.random() > prob
                       else np.random.choice(problem.ALPHABET) for c in X[i, 0]]
                X[i, 0] = "".join(mut)

        return X
```

### Duplicates

Moreover, do not underestimate the importance of filtering out duplicates during the evolution. This helps to make sure a mating produces an offspring that does not already exist in the current population. If it does another mating process is triggered until unique individual do exist.

```python
from pymoo.model.duplicate import ElementwiseDuplicateElimination

class MyDuplicateElimination(ElementwiseDuplicateElimination):

    def is_equal(self, a, b):
        return a.X[0] == b.X[0]

```

### Optimize

Finally, we create an algorithm object with all the modules defined above.

```python
import string
import numpy as np

from pymoo.algorithms.nsga2 import NSGA2
from pymoo.optimize import minimize

algorithm = NSGA2(pop_size=11,
                  sampling=MySampling(),
                  crossover=MyCrossover(),
                  mutation=MyMutation(),
                  eliminate_duplicates=MyDuplicateElimination())

res = minimize(MyProblem(),
               algorithm,
               ('n_gen', 100),
               seed=1,
               verbose=False)

```

```python
from pymoo.visualization.scatter import Scatter
Scatter().add(res.F).show()
```

```python
results = res.X[np.argsort(res.F[:, 0])]
count = [np.sum([e == "a" for e in r]) for r in results[:, 0]]
print(np.column_stack([results, count]))
```

    [['aaaaaaaaaa' 10]
     ['aaaaabaaaa' 9]
     ['abaabaaaaa' 8]
     ['abaaabbaaa' 7]
     ['abbaabaaab' 6]
     ['aabbababab' 5]
     ['bbabaabbba' 4]
     ['ababbabbbb' 3]
     ['bbabbbabbb' 2]
     ['babbbbbbbb' 1]
     ['bbbbbbbbbb' 0]]


text:
```python
# !pip install sympy==1.4
%matplotlib inline
import matplotlib.pyplot as plt
import sympy as sym
from sympy import oo
sym.init_printing()

def fourier_transform(x):
    return sym.transforms._fourier_transform(x, t, w, 1, -1, 'Fourier')
    #return sym.integrate(x*sym.exp(-1j*w*t), (t,0,oo)).evalf()
```

## Definición

La [transformada de Fourier](https://en.wikipedia.org/wiki/Fourier_transform) está definida por

\begin{equation}
X(j \omega) = \int_{-\infty}^{\infty} x(t) \, e^{-j \omega t} \; dt
\end{equation}

donde $X(j \omega) = \mathcal{F} \{ x(t) \}$ se usa como notación de la Transformada de Fourier de la señal $x(t)$. $X(j \omega)$ es el espectro de la señal $x(t)$. El argumento $j \omega$, como exponente de la exponencial, encierra el comportamiento de todas las señales oscilatorias $cos(\omega t)$.

Observe que la forma de la transformada de Fourier corresponde a la forma de la correlación. De esta manera, podría interpretarse como el "parecido" entre la señal $x(t)$ y  $e^{j \omega t}$, es decir entre $x(t)$ y $cos(\omega t)$

La transformada inversa de Fourier $x(t) = \mathcal{F}^{-1} \{ X(j \omega) \}$ se define como

\begin{equation}
x(t) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} X(j \omega) \, e^{j \omega t} \; d\omega
\end{equation}

## Propiedades

### Invertible

\begin{equation}
x(t) = \mathcal{F}^{-1} \left\{ \mathcal{F} \{ x(t) \} \right\}
\end{equation}

Tomando las expresiones de la Transformada de Fourier y la Transformada Inversa de Fourier, se obtiene:

\begin{equation}
\begin{split}
x(t) &= \frac{1}{2 \pi} \int_{-\infty}^{\infty} \underbrace{\int_{-\infty}^{\infty} x(\tau) e^{-j \omega \tau} d\tau}_{X(j \omega)} \; e^{j \omega t} d\omega \\
&= \int_{-\infty}^{\infty} x(\tau) \left( \frac{1}{2 \pi} \int_{-\infty}^{\infty} e^{-j \omega \tau} e^{j \omega t}  d\omega \right) d\tau \\
&= \int_{-\infty}^{\infty} x(\tau) \delta(t - \tau) d\tau = x(t)
\end{split}
\end{equation}

### Linealidad

\begin{equation}
\mathcal{F} \{ A \cdot x_1(t) + B \cdot x_2(t) \} = A \cdot X_1(j \omega) + B \cdot X_2(j \omega)
\end{equation}

Tomando la expresión de la Transformada de Fourier se obtiene:

\begin{equation}
\begin{split}
&= \int_{-\infty}^{\infty} (A \cdot x_1(t) + B \cdot x_2(t)) \, e^{-j \omega t} \; dt \\
&= \int_{-\infty}^{\infty} A \cdot x_1(t) \, e^{-j \omega t} \; dt + \int_{-\infty}^{\infty} B \cdot x_2(t) \, e^{-j \omega t} \; dt \\
&= A \cdot \int_{-\infty}^{\infty} x_1(t) \, e^{-j \omega t} \; dt + B \cdot\int_{-\infty}^{\infty} x_2(t) \, e^{-j \omega t} \; dt
\end{split}
\end{equation}

**Ejemplo - Transformada de Fourier de una señal exponencial causal**

\begin{equation}
x(t) = e^{- \alpha t} \cdot \epsilon(t)
\end{equation}

con $\alpha \in \mathbb{R}^+$

```python
t,w = sym.symbols('t omega', real=True)
a = 4
x = sym.exp(-a * t)*sym.Heaviside(t)
x
```

```python
X = fourier_transform(x)
X
```

```python
plt.rcParams['figure.figsize'] = 5, 2
sym.plot(x, (t,-1,10), ylabel=r'Amp',line_color='blue',legend=True, label = 'x(t)')
sym.plot(sym.re(X), (w,-20,20), ylabel=r'Real',line_color='blue',legend=True, label = 'X(w)')
sym.plot(sym.im(X), (w,-20,20), ylabel=r'Imag',line_color='blue',legend=True, label = 'X(w)')

sym.plot(sym.sqrt(
         (sym.im(X)*sym.im(X)) + (sym.re(X)*sym.re(X))),
         (w,-20,20), ylabel=r'Mag',line_color='blue',legend=True, label = 'X(w)')

```

Observe que:

- $X(\omega)$ es una función definida para todos los valores de $\omega$ y no solamente para los múltiplos enteros de un valor determinado $\omega_0$.

- $X(\omega)$ es una función compleja, es decir que tiene parte imaginaria y parte real. Así, puede expresarse de forma cartesiana ($real + j \cdot imaginario$) o de forma polar ($magnitud \angle ángulo$). El "parecido" entre la señal $x(t)$ con $sin(\omega t)$ se puede apreciar en la magnitud de $X(\omega)$.

- $|X(\omega)|$ tiene un valor máximo en $\omega=0$ y un decaimiento a medida que aumenta $\omega$

Analizando la magnitud de $X(\omega)$

```python
X
```

```python
X_real  = sym.re(X)
X_real
```

```python
X_imag  = sym.im(X)
X_imag
```

```python
X_magn = sym.sqrt(X_real*X_real + X_imag*X_imag).simplify()
X_magn
```

La magnitud de $X(\omega)$ es simétrica respecto a $\omega = 0$. Así, será suficiente analizar solamente un lado del espectro de una señal de tiempo continuo.

**Ejemplo - Transformada de Fourier de una señal exponencial por una senoidal**

\begin{equation}
x(t) = sin(\omega_0 t) \cdot e^{- \alpha t} \cdot \epsilon(t)
\end{equation}

con $\omega \in \mathbb{R}^+$

```python
t,w = sym.symbols('t omega', real=True)
w0 = 10
x1 = sym.sin(w0 * t)*sym.exp(-2*t)*sym.Heaviside(t)
x1
```

```python
X1 = fourier_transform(x1)
X1
```

```python
plt.rcParams['figure.figsize'] = 5, 2
sym.plot(x1, (t,-2,5), ylabel=r'Amp',line_color='blue',legend=True, label = 'x1(t)')
#sym.plot(sym.re(X), (w,-20,20), ylabel=r'Real',line_color='blue',legend=True, label = 'X(w)')
#sym.plot(sym.im(X), (w,-20,20), ylabel=r'Imag',line_color='blue',legend=True, label = 'X(w)')

sym.plot(sym.sqrt(
         (sym.im(X1)*sym.im(X1)) + (sym.re(X1)*sym.re(X1))),
         (w,-60,60), ylabel=r'Mag',line_color='blue',legend=False, label = 'X1(w)')

```

Observe que:
- $x1(t)$ corresponde a $x(t)$ multiplicada con una función senoidal de frecuencia angular de $10$ rad/seg.

- Al igual que con la magnitud de $X(\omega)$, la magnitud de $X1(\omega)$ decae con $\omega$, sin embargo, hay un pico en $\omega=10$ que se relaciona con la senoidal.

**Ejemplo - Transformada de Fourier de senoidales y combinaciones causales**

\begin{equation}
x(t) = sin(\omega_0 t) \cdot \epsilon(t)
\end{equation}

con $\omega \in \mathbb{R}^+$

```python
w1 = 10
w2 = 5
x2_1 = sym.sin(w1 * t)*sym.Heaviside(t)

x2_2 = sym.sin(w2 * t)*sym.Heaviside(t)

x2 = x2_1 + x2_2
x2
```

```python
X2_1 = fourier_transform(x2_1)
X2_1

```

```python
X2_2 = fourier_transform(x2_2)
X2_2
```

```python
(X2_1+X2_2).simplify()
```

```python
X2 = fourier_transform(x2)
X2
```

```python
plt.rcParams['figure.figsize'] = 5, 2
gt2_1 = sym.plot(x2_1, (t,-1,5), ylabel=r'Amp',line_color='blue',legend=True, label = 'x2_1(t)', show = False)
gt2_2 = sym.plot(x2_2, (t,-1,5), ylabel=r'Amp',line_color='green',legend=True, label = 'x2_2(t)', show = False)
gt2 = sym.plot(x2, (t,-1,5), ylabel=r'Amp',line_color='red',legend=True, label = 'x2(t)', show = False)
gt2.extend(gt2_1)
gt2.extend(gt2_2)
gt2.show()

```

```python

plt.rcParams['figure.figsize'] = 6, 3

gw2_1 = sym.plot(sym.sqrt(
         (sym.im(X2_1)*sym.im(X2_1)) + (sym.re(X2_1)*sym.re(X2_1))),
         (w,0,14), ylabel=r'Mag',line_color='blue',legend=False, label = 'X2_1(w)',show = False)
gw2_2 = sym.plot(sym.sqrt(
         (sym.im(X2_2)*sym.im(X2_2)) + (sym.re(X2_2)*sym.re(X2_2))),
         (w,0,14), ylabel=r'Mag',line_color='green',legend=False, label = 'X2_2(w)',show = False)
gw2 = sym.plot(sym.sqrt(
         (sym.im(X2)*sym.im(X2)) + (sym.re(X2)*sym.re(X2))),
         (w,0,14), ylabel=r'Mag',line_color='red',legend=False, label = 'X2(w)',show = False)

gw2.extend(gw2_1)
gw2.extend(gw2_2)
gw2.show()
```

```python
X2
```

En la gráfica anterior se observa el efecto de la superposición lineal de los espectros de las dos señales senoidales.

**Ejercicio**

Analice el espectro de

\begin{equation}
x(t) = (sin(\omega_0 t) + e^{-2t}) \cdot \epsilon(t)
\end{equation}

con $\omega_0 \in \mathbb{R}^+$

```python
t,w = sym.symbols('t omega', real=True)
w0 = 15
a = 2
x4 = (sym.sin(w0 * t) + sym.exp(-a*t))*sym.Heaviside(t)
x4
```

```python

```

### Dualidad

Observe que la **Transformada de Fourier** y la **Transformada Inversa de Fourier** tienen formas parecidas.

\begin{align}
X(\omega) &= \int_{-\infty}^{\infty} x(t) \, e^{-j \omega t} \; dt \\
x(t) &= \frac{1}{2 \pi} \int_{-\infty}^{\infty} X(j \omega) \, e^{j \omega t} \; d\omega
\end{align}

La principal diferencia está en el factor de normalización $2 \pi$ y el signo de la exponencial.

Suponga que:

\begin{equation}
x_2(\omega) = \mathcal{F} \{ x_1(t) \}
\end{equation}

Puede pensarse que:

\begin{equation}
x_2(t) = x_2(\omega) \big\vert_{\omega=t}
\end{equation}

Entonces

\begin{equation}
\mathcal{F} \{ x_2(t) \} = \int_{-\infty}^{\infty} x_2(\omega) \big\vert_{\omega=t} \, e^{-j \omega t} \; dt
\end{equation}

Esta tiene la forma de **Transformada de Fourier**, pero la función interna tiene a $\omega$ como variable, esto indica que la integral se trata como una **Transformada inversa**.  Así, para volver a la **Transformada de Fourier**, se debe multiplicar por $2\pi$ y quitar el signo de la exponencial del kernel de transformación.

\begin{equation}
\mathcal{F} \{ x_2(t) \} = 2 \pi \cdot x_1(- \omega)
\end{equation}

Esta propiedad permite llevar los análisis desde el dominio de las frecuencias al dominio del tiempo y viceversa.

## Teoremas

Retomando la transformada

\begin{equation}
X(j \omega) = \int_{-\infty}^{\infty} x(t) \, e^{-j \omega t} \; dt
\end{equation}

### Derivadas

Dadas una señal $x(t)$ y su derivada respecto al tiempo $\frac{d x(t)}{dt}$, y conocida la **Transformada de Fourier** $X(\omega)$ :

\begin{equation}
\mathcal{F} \left\{ \frac{d x(t)}{dt} \right\}  = \int_{-\infty}^{\infty} \frac{d x(t)}{dt} \, e^{-j \omega t} \; dt
\end{equation}

La integral se puede resolver por partes:

\begin{equation}
\begin{split}
\mathcal{F} \left\{ \frac{d x(t)}{dt} \right\} &= x(t) \cdot e^{-j \omega t} \big\vert_{-\infty}^{\infty} - \int_{-\infty}^{\infty} x(t) (-j \omega) e^{-j \omega t} \; dt \\
&= j \omega \int_{-\infty}^{\infty} x(t) e^{-j \omega t} \; dt \\
&= j \omega X(\omega) \
\end{split}
\end{equation}

\begin{equation}
\frac{d x(t)}{dt} = \frac{d \delta(t)}{dt} * x(t)
\end{equation}

La principal aplicación está en la transformación de ecuaciones diferenciales.

** Ejemplo**

\begin{equation}
2y(t) + 2 \frac{dy}{dt} - x(t) = 0
\end{equation}

Aplicando la **Transformada de Fourier** con sus propiedades se obtiene:
\begin{equation}
2Y(\omega) + 2 j \omega Y(\omega) - X(\omega) = 0
\end{equation}

Observe que en el modelo en el dominio del tiempo (ecuación diferencial) no es posible despejar una expresión equivalente a $\frac{x(t)}{y(t)}$. Por su parte, usando el modelo en el dominio de la frecuencia, se obtiene:

\begin{equation}
Y(\omega)(2 + 2 j \omega ) = X(\omega)
\end{equation}

\begin{equation}
\frac{Y(\omega)}{X(\omega)} = \frac{1}{2+2j\omega} = F(\omega)
\end{equation}

Esta relación es conocida como **Función de transferencia** y representa el efecto que tiene el sistema sobre una señal de entrada que en el caso de la transformada de Fourier es senoidal.

```python
# La función de transferencia
F = 1 / (2 +1j*2*w  )

F
```

```python
plt.rcParams['figure.figsize'] = 5, 2
sym.plot(sym.Abs(F), (w,0,10), ylabel=r'Mag',line_color='blue',legend=True, label = 'F(w)', show = True)

```

La respuesta de un sistema ante una entrada senoidal de frecuencia específica $\omega$ se determina por el valor complejo que toma $F(\omega)$. Por ejemplo, en la frecuencia $\omega = 1$, $F(1) = \frac{1}{2+2j}$

```python
F1 = F.subs(w,1)
F1

```

```python
magF1 = sym.Abs(F1)
magF1
```

Así, si el sistema es excitado con un seno de frecuencia 1, con amplitud 1, la salida será un seno de amplitud 0.35 y tendrá un desfase definido por el ángulo de la función de transferencia.

\begin{equation}
0.35 sin(1t + ang)
\end{equation}

### Ejercicio

Deduzca cómo se calcula el ángulo de fase que introduce el sistema.

```python
F
```

```python
F.subs(w,3)
```

```python
sym.re(F)
```

```python
sym.im(F)
```

```python
sym.sqrt(sym.re(F)**2 +
     sym.im(F)**2  )
```

```python
sym.atan(sym.im(F)/sym.re(F))
```

```python
sym.Abs(F)
```

```python
sym.arg(F)
```

```python
plt.rcParams['figure.figsize'] = 5, 2
sym.plot(sym.Abs(F), (w,0,10), ylabel=r'Mag',line_color='blue',legend=True, label = 'F(w)', show = True)
sym.plot(sym.arg(F), (w,0,10), ylabel=r'arg',line_color='blue',legend=True, label = 'F(w)', show = True)

```

```python

```


text:
### DEMDP01
# Timber Harvesting Model - Cubic Spline Approximation

Profit maximizing owner of a commercial tree stand must decide when to clearcut the stand.

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from compecon import DPmodel, BasisSpline, demo, NLP
```

### Model Parameters
Assuming that the unit price of biomass is $p=1$, the cost to clearcut-replant is $\kappa=0.2$, the stand carrying capacity $s_{\max}=0.5$, biomass growth factor is $\gamma=10\%$ per period, and the annual discount factor $\delta=0.9$.

```python
price = 1.0
kappa = 0.2
smax  = 0.5
gamma = 0.1
delta = 0.9
```

### State Space
The state variable is the stand biomass, $s\in [0,s_{\max}]$.

Here, we represent it with a cubic spline basis, with $n=200$ nodes.

```python
n = 200
basis = BasisSpline(n, 0, smax, labels=['biomass'])
```

### Action Space
The action variable is $j\in\{0:\text{'keep'},\quad 1:\text{'clear cut'}\}$.

```python
options = ['keep', 'clear-cut']
```

### Reward function
If the farmer clears the stand, the profit is the value of selling the biomass $ps$ minus the cost of clearing and replanting $\kappa$, otherwise the profit is zero.

```python
def reward(s, x, i , j):
    return (price * s - kappa) * j
```

### State Transition Function
If the farmer clears the stand, it begins next period with $\gamma s_{\max}$ units. If he keeps the stand, then it grows to $s + \gamma (s_{\max} - s)$.

```python
def transition(s, x, i, j, in_, e):
    if j:
        return np.full_like(s, gamma * smax)
    else:
        return s + gamma * (smax - s)
```

### Model Structure
The value of the stand, given that it contains $s$ units of biomass at the beginning of the period, satisfies the Bellman equation

\begin{equation} V(s) = \max\left\{(ps-\kappa) + \delta V(\gamma s_{\max}),\quad \delta V[s+\gamma(s_{\max}-s)]  \right\}   \end{equation}

To solve and simulate this model, use the CompEcon class ```DPmodel```.

```python
model = DPmodel(basis, reward, transition,
                discount=delta,
                j=options)

S = model.solve()
```

    Solving infinite-horizon model collocation equation by Newton's method
    iter change       time
    ------------------------------
       0       3.1e-01    0.0000
       1       1.1e-01    0.0186
       2       1.1e-02    0.0221
       3       4.6e-03    0.0221
       4       1.9e-16    0.0378
    Elapsed Time =    0.04 Seconds

The ```solve``` method retuns a pandas ```DataFrame```, which can easily be used to make plots. To see the first 10 entries of `S`, we use the `head` method

```python
S.head()
```

      biomass
      value
      resid
      j*
      value[keep]
      value[clear-cut]

      biomass

      0.000000
      0.000000
      0.067287
      1.387779e-17
      keep
      0.067287
      -0.132713

      0.000250
      0.000250
      0.067320
      -7.491474e-09
      keep
      0.067320
      -0.132462

      0.000500
      0.000500
      0.067353
      -5.953236e-09
      keep
      0.067353
      -0.132212

      0.000750
      0.000750
      0.067387
      -1.332984e-09
      keep
      0.067387
      -0.131962

      0.001001
      0.001001
      0.067420
      7.309628e-10
      keep
      0.067420
      -0.131712

## Analysis

To find the biomass level at which it is indifferent to keep or to clear cut, we interpolate as follows:

```python
scrit = NLP(lambda s: model.Value_j(s).dot([1,-1])).broyden(0.0)[0]
vcrit = model.Value(scrit)
print(f'When the biomass is {scrit:.5f} its value is {vcrit:.5f} regardless of the action taken by the manager')
```

    When the biomass is 0.30669 its value is 0.17402 regardless of the action taken by the manager

### Plot the Value Function and Optimal Action

```python
fig1 = plt.figure(figsize=[8,5])
gs = gridspec.GridSpec(2, 1,height_ratios=[3.5, 1])

S[['value[keep]', 'value[clear-cut]']].plot(ax= plt.subplot(gs[0]))
plt.title('Action-Contingent Value Functions')
plt.xlabel('')
plt.ylabel('Value of Stand')
plt.xticks([])

ymin,ymax = plt.ylim()
plt.vlines(scrit,ymin, vcrit,linestyle=':')
demo.annotate(scrit, ymin,'$s^*$',ms=0)

S['j*'].cat.codes.plot(ax=plt.subplot(gs[1]))
plt.title('Optimal Action')
plt.ylabel('Action')
plt.ylim([-0.25,1.25])
plt.yticks([0,1],options)
plt.legend([])
```

### Plot Residuals

```python
S['resid2'] = 100*S.resid / S.value

fig2 = demo.figure('Bellman Equation Residual','','Percent Residual')
plt.plot(S['resid2'])
plt.hlines(0,0,smax,'k')
```

###  Simulation

The path followed by the biomass is computed by the ```simulate()``` method. Here we simulate 32 periods starting with a biomass level $s_0 = 0$.

```python
H = model.simulate(32, 0.0)

fig3 = demo.figure('Timber harvesting simulation','Period','Biomass')
plt.plot(H['biomass'])

```

```python
#demo.savefig([fig1, fig2, fig3])
```

max_stars_repo_path: notebooks/dp/01 Timber Harvesting Model - Cubic Spline Approximation.ipynb

max_stars_repo_name: daniel-schaefer/CompEcon-python

max_issues_repo_name: daniel-schaefer/CompEcon-python

max_forks_repo_name: daniel-schaefer/CompEcon-python