policy_gradient_proof

Proof of the Policy Gradient Theorem

一些说在证明前面的话

自己看了Sutton&Barto的证明(Sutton书中的13.2节的一个box的证明过程)，参考了其前四步的推理，觉得是非常清晰且正确的，但是到第五步将下面推导的第四步标红的$\nabla v_{\pi }$再次展开并写成这种格式时觉得非常难理解

$$\nabla v_{\pi }(s)=\sum _{x\in \mathcal{S}}\left(\sum _{k=0}^{\infty }\mathrm{Pr}(s\to x,k,\pi )\right)\sum _a\nabla \pi (a|s)q_{\pi }(s,a)$$

• 这里的$\mathrm{Pr}(s\to x,k,\pi )$是一个概率，表示从状态$s$经过$k$步转移到状态$x$的概率，和之前的证明average reward中的$p_{\pi }^{(k)}(x|s)$其实是一样的
• 但是从接下来我们所看到的第四步的推导中，我们发现并没有办法直接推导出上述形式，在这里笔者决定第五步之后采取了另一种做法——矩阵形式的推导，这样可以更加清晰地展示出推导的过程，最终会证明上述这个式子是对的（从级数展开的角度去理解）

Step 1: Unroll $\nabla v_{\pi }(s)$ to iterative form(Refer Sutton&Barto)

We can prove the policy gradient theorem from first principles.

• $\pi$ is a function of $\theta$, means that we use $\pi (a|s,\theta )\to \pi (a|s)$ for simplicity
• all gradients are with respect to $\theta$, means that we use ${\nabla }_{\theta }\to \nabla$ for simplicity

$$\begin{array}{rlrl}\nabla v_{\pi }(s)& =\nabla \left[\sum _a\pi (a|s)q_{\pi }(s,a)\right],\forall s\in \mathcal{S}& & \\ & =\sum _a\left[\nabla \pi (a|s)q_{\pi }(s,a)+\pi (a|s)\nabla q_{\pi }(s,a)\right]& & \text{(product rule of calculus)}\\ & =\sum _a\left[\nabla \pi (a|s)q_{\pi }(s,a)+\pi (a|s)\nabla \sum _{s^{\prime },r}p(s^{\prime },r|s,a)[r+v_{\pi }(s^{\prime })]\right]& & (\text{Eq. 1})\\ & =\sum _a\left[\nabla \pi (a|s)q_{\pi }(s,a)+\pi (a|s)\sum _{s^{\prime }}p(s^{\prime }|s,a)\nabla v_{\pi }(s^{\prime })\right]& & (\text{Eq. 2 and r is a constant})\end{array}$$

where

$$\begin{array}{rlrl}\text{Eq. 1:}& & q_{\pi }(s,a)& =\sum _{s^{\prime },r}p(s^{\prime },r|s,a)[r+v_{\pi }(s^{\prime })]\\ \text{Eq. 2:}& & p(s^{\prime }|s,a)& =\sum _{r\in \mathcal{R}}p(s^{\prime },r|s,a)\end{array}$$

Step 2: Turn $\nabla v_{\pi }(s)$ into matrix form

Refer to our former notes on the probability&matrix, we know that any conditional probability can turn into a matrix easily, just do it:

1. $\nabla v_{\pi }(s)\to \nabla v_{\pi }$ is a (n,1) vector
   • the shape is same as $v_{\pi }$
   • where $n$ is the number of states
2. $\nabla \pi (a|s)\to \nabla \pi$ is a (n,m) matrix
   • the shape is same as $\pi$
   • where $m$ is the number of actions
3. $q_{\pi }(s,a)\to q_{\pi }$ is a (n,m) matrix
   • the shape is same as $\pi$
4. $p(s^{\prime }|s,a)\to P_s$ is a (n,m,n) tensor

Now, we can rewrite the above last equation in matrix form.

$$\nabla v_{\pi }=\nabla \pi @q_{\pi }+\pi @P_s\nabla v_{\pi }$$

where $@$ denotes batch matrix multiplication(bmm).

We demonstrate the above matrix equation in a more detailed form.

• first, we caculate the $\nabla \pi @q_{\pi }$ → (n,1)
  • $\nabla \pi$ is a (n,m) matrix, we unsequeeze it to (n,1,m) tensor
  • $q_{\pi }$ is a (n,m) matrix, we unsequeeze it to (n,m,1) tensor
  • then, we do bmm, the result is a (n,1,1) tensor
  • finally, we squeeze it to (n,1) vector

Why we do bmm $@$ here?

$$\sum _a\nabla \pi (a|s)q_{\pi }(s,a)\Rightarrow \nabla \pi @q_{\pi }$$

We notice ${\sum }_a$ is the dot-product-sum of the two tensors $\nabla \pi$ (n,m) dim=1 and $q_{\pi }$ (n,m) dim=1, and we should also return a tensor shape like $v_{\pi }$, which is (n,1), so we get (n,1,m)@(n,m,1) = (n,1,1), then we squeeze it to (n,1) vector

• second, we caculate the $\pi @P_s\nabla v_{\pi }$ → (n,1)
  • $\pi$ is a (n,m) matrix, we unsequeeze it to (n,1,m) tensor
  • $P_s$ is a (n,m,n) tensor, we just use it directly
  • $\nabla v_{\pi }$ is a (n,1) vector, use it directly(broadcast)
  • then, we do bmm, $\pi @P_s$ is a (n,1,n) tensor, squeeze it to (n,n)
  • finally, just do matrix multiplication, the result is a (n,1) vector

Why we do bmm $@$ here?

$$\sum _a\pi (a|s)\sum _{s^{\prime }}p(s^{\prime }|s,a)\nabla v_{\pi }(s^{\prime })\Rightarrow \pi @P_s\nabla v_{\pi }$$

• ${\sum }_{s^{\prime }}$ is the dot-product-sum of $P_s$ (n,m,n) dim=2 and $\nabla v_{\pi }$ (n,1) dim=0, actually, we don’t need to do bmm here, just matmul and we get (n,m,n)x(n,1) = (n,m,1), then we squeeze it to (n,m) tensor
• ${\sum }_a$ is the dot-product-sum of the two tensors $\pi$ (n,m) dim=1 and the returned tensor (n,m) from the last step, so we get (n,1,m)@(n,m,1) = (n,1,1), then we squeeze it to (n,1) vector
• Note the order of ${\sum }_{s^{\prime }}$ and ${\sum }_a$ is not important, because the result is the same

Step 3: Solve the matrix equation

Directly unroll the matrix equation

$$\begin{array}{rlrl}\nabla v_{\pi }& =\nabla \pi @q_{\pi }+\pi @P_s\nabla v_{\pi }& & \\ & =\nabla \pi @q_{\pi }+\pi @P_s(\nabla \pi @q_{\pi }+\pi @P_s\nabla v_{\pi })& & (\text{substitute ∇vπ})\\ & =(I+\pi @P_s)\nabla \pi @q_{\pi }+(\pi @P_s{)}^2\nabla v_{\pi }& & \\ & & & \\ & =(I+P_{\pi })\nabla \pi @q_{\pi }+P_{\pi }^2\nabla v_{\pi }& & (P_{\pi }=\pi @P_s)\\ & =(I+P_{\pi }+\cdots +P_{\pi }^{k-1})\nabla \pi @q_{\pi }+P_{\pi }^k\nabla v_{\pi }& & (\text{expand to k steps})\\ & =(I-P_{\pi }^k)(I-P_{\pi }{)}^{-1}\nabla \pi @q_{\pi }+P_{\pi }^k\nabla v_{\pi }& & (\text{use geometric series})\end{array}$$

什么是geometric series(几何级数)?

就是指等比数列前$n$项和，这里写成了矩阵形式，其实就是

$$\begin{array}{rl}{S}_k& =I+P_{\pi }+\cdots +P_{\pi }^{k-1}\\ \Rightarrow S_{k+1}& =S_k+P_{\pi }^k=I+S_k{P}_{\pi }\\ \Rightarrow S_k& =(I-P_{\pi }^k)(I-P_{\pi }{)}^{-1}\end{array}$$

We notice that from Appendix

$$\begin{array}{rl}{d}_{\pi }^T\underset{k\to \infty }{\lim }{P}_{\pi }^k& =d_{\pi }^T{P}_{\pi }=d_{\pi }^T\\ \Rightarrow \underset{k\to \infty }{\lim }{P}_{\pi }^k& =C\end{array}$$

Then we have another beautiful iterative form here:

$$\nabla v_{\pi }=(I-C)(I-P_{\pi }{)}^{-1}\nabla \pi @q_{\pi }+C\nabla v_{\pi }$$

And we can also obtain the final matrix form from above:

$$\begin{array}{rl}(I-C)\nabla v_{\pi }& =(I-C)(I-P_{\pi }{)}^{-1}\nabla \pi @q_{\pi }\\ \Rightarrow v_{\pi }& =(I-P_{\pi }{)}^{-1}\nabla \pi @q_{\pi }\end{array}$$

理解最开始 说在证明前面的话中的式子

上面这个式子其实就是最开头的式子的matrix form，其中，值得注意的是（也是比较难以一眼看出的就是下面这个概率到matrix的转换）

$$\sum _{k=0}^{\infty }\mathrm{Pr}(s\to x,k,\pi )\Rightarrow (I-P_{\pi }{)}^{-1}$$

从级数展开的角度去理解就很清晰了：

$$(I-P_{\pi }{)}^{-1}=I+P_{\pi }+P_{\pi }^2+\cdots$$

Step 4: Obtain the final form(Refer Sutton&Barto)

Finally, we can obtain the final form of the policy gradient theorem.

• matrix form $$\nabla v_{\pi }=(I-P_{\pi }{)}^{-1}\nabla \pi @q_{\pi }$$
• element-wise form $$\nabla v_{\pi }(s)=\sum _{x\in \mathcal{S}}\left(\sum _{k=0}^{\infty }\mathrm{Pr}(s\to x,k,\pi )\right)\sum _a\nabla \pi (a|x)q_{\pi }(x,a)$$

It is then immediate that

$$\begin{array}{rlrl}J(\theta )& \doteq \nabla v_{\pi }(s_0)& & \\ & =\sum _s\left(\sum _{k=0}^{\infty }\mathrm{Pr}(s_0\to s,k,\pi )\right)\sum _a\nabla \pi (a|s)q_{\pi }(s,a)& & \\ & =\sum _s\eta (s)\sum _a\nabla \pi (a|s)q_{\pi }(s,a)& & (\text{time steps spent})\\ & =\sum _s\sum _{s^{\prime }}\frac{\eta (s^{\prime })}{{\sum }_{s^{\prime }}\eta (s^{\prime })}\eta (s)\sum _a\nabla \pi (a|s)q_{\pi }(s,a)& & (\sum _{s^{\prime }}\frac{\eta (s^{\prime })}{{\sum }_{s^{\prime }}\eta (s^{\prime })}=1)\\ & =\sum _{s^{\prime }}\eta (s^{\prime })\sum _s\frac{\eta (s)}{{\sum }_{s^{\prime }}\eta (s^{\prime })}\sum _a\nabla \pi (a|s)q_{\pi }(s,a)& & (\text{rearrange the sum})\\ & =\sum _{s^{\prime }}\eta (s^{\prime })\sum _s\mu (s)\sum _a\nabla \pi (a|s)q_{\pi }(s,a)& & (\text{on-policy distribution})\\ & \propto \sum _s\mu (s)\sum _a\nabla \pi (a|s)q_{\pi }(s,a)& & (\text{Q.E.D.})\end{array}$$

where

• $\eta (s)$ is the time steps spent(or number of visits) in state $s$, details see Appendix
• $\mu (s)$ is just the sum-to-one normalized $\eta (s)$, we call it on-policy distribution
• note that the $s^{\prime }$ here does not mean the next state, it is just a variable notation to distinguish the sum
• ${\sum }_{s^{\prime }}\eta (s^{\prime })=C$, just a constant, so we can ignore it

What is Q.E.D?

有几种叫法：

1. 希腊语：ὅπερ ἔδει δεῖξαι（hoper edei deixai）欧几里得和阿基米德用过
2. 英语：Q.E.D (quod erat demonstrandum) 希腊语的翻译，Sutton&Barto用过
3. 中文：证毕，我用过😏

TL;DR

Using the matrix form, we can quickly summarize the above proof process. (Note that the $\nabla$ here is just the same as the operator $\mathrm{d}$, does not represent the derivative of $\theta$, in short $\nabla \ne {\nabla }_{\theta }$)

$$\begin{array}{rlrl}\nabla v_{\pi }& =\nabla (\pi @q_{\pi })& & (v_{\pi }=\pi @q_{\pi })\\ & =\nabla \pi @q_{\pi }+\pi @\nabla q_{\pi }& & (\text{product rule})\\ & =\nabla \pi @q_{\pi }+\gamma \pi @P_s\nabla v_{\pi }& & (q_{\pi }=P_{r}R+\gamma P_s@v_{\pi })\\ & =\nabla \pi @q_{\pi }+\gamma P_{\pi }\nabla v_{\pi }& & (P_{\pi }=\pi @P_s)\\ & =(I-\gamma P_{\pi }{)}^{-1}\nabla \pi @q_{\pi }& & (\text{rearrange})\\ & =\left(\sum _{k=0}^{\infty }(\gamma P_{\pi }{)}^k\right)\nabla \pi @q_{\pi }& & (\text{unroll})\end{array}$$

Conclusion

From the Appendix, we notice that

• without-discounting $\gamma =1$ case: each row of the normalized $\left({\sum }_{k=0}^{\infty }{P}_{\pi }^k\right)$ is same, in other word, each row is the stationary distribution $d_{\pi }$. Therefore, we can just take one value from the vector $\nabla v_{\pi }$ to represent the scalar objective function $J(\theta )$.

  $$\begin{array}{rl}\nabla J(\theta )& =\nabla v_{\pi }(s)\\ & =k{d}_{\pi }^T\nabla \pi @q_{\pi }\\ & \propto d_{\pi }^T\nabla \pi @q_{\pi }\end{array}$$
  • $d_{\pi }$ is the stationary distribution, it is a (n,1) vector
  • $\nabla \pi @q_{\pi }$ is a (n,1) vector as we mentioned before
  • $k$ is just a constant, it will be very large because $k\to \infty$

• discounting $\gamma <1$ case: each row of the normalized $\left({\sum }_{k=0}^{\infty }(\gamma P_{\pi }{)}^k\right)$ is different, so we can not just take one value from the vector $\nabla v_{\pi }$ to represent the scalar objective function $J(\theta )$. We can use a $\nabla {\bar{v}}_{\pi }$ to represent the scalar objective function $J(\theta )$.

  $$\begin{array}{rl}\nabla J(\theta )& =\nabla d^T{v}_{\pi }(s)\\ & =d^T\mathrm{M}C\nabla \pi @q_{\pi }\\ & \propto d^T\mathrm{M}\nabla \pi @q_{\pi }\end{array}$$
  • $d$ is the current distribution of state, it is a (n,1) vector
  • each row of $\mathrm{M}$ is a on-policy distribution, it is a (n,n) matrix
  • $C$ is just a constant (n,1) vector, each row is same, we can take it as a scalar $C=1/(1-\gamma )$

• unified form

  $$\boxed{\nabla J(\theta )\propto \mu \nabla \pi @q_{\pi }}$$
  • $\mu =d^T\mathrm{M}$ is the mixed on-policy distribution, it is a (1,n) vector
  • $\gamma =1$, each row of $\mathrm{M}$ is $d_{\pi }$, so $\mu =d_{\pi }^T$
  • $\gamma <1$, each row of $\mathrm{M}$ is different, so $\mu =d^T\mathrm{M}$

一些思考

自己一开始没留意着原来discount和没有discount的区别会导致最后的结果的差异那么大，确实是自己忽略这部分原因思考了，下面Appendix中的内容说明了这个常数C的值其实，在下一节做优化（梯度上升）时自己觉得可以稍微考虑一下C的因素在里面，当$\gamma$非常接近1时就考虑用比例来做；当$\gamma$小于1时，可以考虑用原来的式子来做，即把$C=1/(1-\gamma )$代入到$\nabla J(\theta )$中，在梯度上升时可以考虑这个常数，就不用啥都归入$\alpha$中来调参比较困难了。也即我们之后可以考虑的正确梯度形式为（注意这里是等号，而不是成比例）：

$$\nabla J(\theta )=C\mu \nabla \pi @q_{\pi }$$

而且这个结果还可以从另外一个地方得到互相验证：

1. 上面通篇的推导我们都是用的metric ${\bar{v}}_{\pi }$, 而在之前的章节中我们证明了另一个metric ${\bar{r}}_{\pi }=(1-\gamma ){\bar{v}}_{\pi }$，也就是说如果用这个${\bar{r}}_{\pi }$作为metric，那么我们的policy gradient就变成了$\nabla J(\theta )=\mu \nabla \pi @q_{\pi }$（还是那句话这里是等号）；
2. 如果假设我们很笨，一开始没有推出第二个metric，我们可以仍然只用metric ${\bar{v}}_{\pi }$进行思考，注意力放在$C=1/(1-\gamma )$上，$\gamma =0$所导致的结果就是我们只着眼于眼前的奖励（回想刚刚引入$\gamma$时的定义discounted factor），注意到$\gamma =0$的metric ${\bar{v}}_{\pi }$和metric ${\bar{r}}_{\pi }$其实是一个式子的policy gradient

另外自己觉得Sutton书中13.3节的REINFORCE伪代码中在最后的梯度部分引入了${\gamma }^t$这个是有问题的，这样的副作用导致Figure 13.1中的曲线中的$\alpha =2^{-13}$，这个值过于小了。在Exercise 13.2的答案中也提及了这里unclear，自己的猜测：由于没考虑常数C的影响导致了学习率需要设置得比较大($\gamma =0.99$的话，本来的alpha=0.1就得设成alpha=10才会有效果了，但一般大家都默认alpha<1)，所以他们才引入了${\gamma }^t$这个因素 ⇒ 正确做法应该是考虑C, 之后实践一下看看效果

Appendix

Notations about long-run behavior

We define some notations about the long-run($k\to \infty$) behavior for the above proof.

$$\begin{array}{rlrl}\mathrm{H}& \doteq \sum _{k=0}^{\infty }(\gamma P_{\pi }{)}^k& & (\text{the upper η})\\ C& \doteq \sum _{s^{\prime }}\mathrm{H}(s^{\prime })=\sum _{\dim =1}\mathrm{H}& & (\text{the constant})\\ \mathrm{M}& \doteq \mathrm{H}/C& & (\text{the upper μ})\end{array}$$

• $\mathrm{H}$ is a (n,n) matrix, $\eta$ is a (1,n) vector
  • $\mathrm{H}(s,s^{\prime })$ is time steps spent (number of visits) starting from $s$ and ending at $s^{\prime }$
  • $\eta (s^{\prime })$ is time steps spent (number of visits) (starting from default $s_0$) ending at $s^{\prime }$
• $\mathrm{M}$ is a (n,n) matrix, $\mu$ is a (1,n) vector
  • $\mathrm{M}(s,s^{\prime })$ is the probability of visits starting from $s$ and ending at $s^{\prime }$
  • $\mu (s^{\prime })$ is the probability of visits (starting from default $s_0$) ending at $s^{\prime }$
• $C$ is a (n,1) vector
  • $C(s)$ is the total visits starting from $s$
  • We find that
    • $\gamma =1$: $C(s_1)=C(s_2)=\cdots =C(s_n)=k-1={\sum }_{i=1}^{k-1}1$
      • $\boxed{\underset{k\to \infty }{\lim }C(s)=\infty }$
    • $\gamma <1$: $C(s_1)=C(s_2)=\cdots =C(s_n)=(1-{\gamma }^{k-1})/(1-\gamma )={\sum }_{i=1}^{k-1}{\gamma }^{i-1}$
      • $\boxed{\underset{k\to \infty }{\lim }C(s)=1/(1-\gamma )}$
  • where $k$ is just the number of iterations or total visits(usaually a large number because $k\to \infty$)

The on-policy distribution(Refer Sutton&Barto)

Let

• $d(s)$ denote the probability that an episode begins in each state $s$
• $\eta (s)$ denote the number of time step spent, on average, in state $s$ in a single episode

Time is spent in a state $s$,

• if episodes start in $s$
• or if transitions are made into $s$ from a preceding(先前的) $\bar{s}$

Then we have the following equations which can be solved for the expected number of visits $\eta (s)$

$$\eta (s)=d(s)+\gamma \sum _s\eta (\bar{s})\sum _a\pi (a|\bar{s})p(s|\bar{s},a),\forall s\in \mathcal{S}$$

• without discounting case: $\gamma =1$
• discounting case: $\gamma <1$

The on-policy distribution is

$$\mu (s)=\frac{\eta (s)}{{\sum }_{s^{\prime }}\eta (s^{\prime })},\forall s\in \mathcal{S}$$

Matrix form:

• Iterative form: $$\eta =d+\gamma \eta \pi @P_s=d+\gamma \eta P_{\pi }$$
• Solution form: $$\eta =(I-\gamma P_{\pi }{)}^{-1}d$$
• And $\mu =\eta /\sum \eta$
• where $\eta ,d,\mu$ is a (n,1) vector

一些理解

这里的主要思想其实就是expected number of visits, $\eta (s)$表示的是state $s$被访问次数的期望，计算他的方法也很简单，是从递归的角度去想的，我们现在把$\eta (s)$理解成是 probability x 1(time-step) spent in a state $s$, 那么我们要做的就是简单的将所有概率加起来就行了，也就是两部分组成

1. 当前state的概率d(s)x1
2. 所有可以跳到当前的state的其他state的被访问次数 x 状态转移概率 x 1并求和（因为假设$\bar{s}$被访问k次，$P_{\pi }(s|\bar{s})=0.1$，那么当然s也会有这部分0.1k次被访问到才对）

Convergence of the long-run behavior

Just modify the example in Appendix, and we will get the same result d_pi (but matrix form)

The following code demonstrates that the three is same to represent the stationary distribution.

• $P_{\pi }^k$ is the transition matrix multiplied $k$ times
• $(I-P_{\pi }{)}^{-1}$ is the inverse of the transition matrix, we normalize it by the sum of each row
• $I+P_{\pi }+P_{\pi }^2+\cdots$ is the geometric series of the transition matrix, we normalize it by the sum of each row

The three results are the same, but

• the first one is the most efficient, easy to calculate
• the second one is the most accurate, but hard to calculate
• the third one is more intuitive, but need a lot of iterations

Talk is cheap, show the code.

# File: P_pi_k.py
import numpy as np
 
def method_multiply(P_pi: np.ndarray, k: int):
    res = np.eye(P_pi.shape[0])
    for i in range(k):
        res = res @ P_pi
    return res
 
def inverse_norm_direct(M: np.ndarray):
    inv = np.linalg.inv(M)
    print(np.sum(inv, axis=1))
    return inv / np.sum(inv, axis=1)
 
def inverse_norm_iterative(M: np.ndarray):
    I = np.eye(M.shape[0])
    inv = np.zeros(M.shape)
    cnt = 0
    # iterate until convergence
    while True:
        cnt += 1
        inv_new = I + inv @ (I - M)
        inv_new = inv_new
        if np.allclose(inv, inv_new):
            break
        inv = inv_new
    print(f"converged after {cnt} iterations")
    print(np.sum(inv, axis=1))
    return inv / np.sum(inv, axis=1)
 
 
P_pi = np.array([
    [0.3, 0.1, 0.6, 0],
    [0.1, 0.3, 0, 0.6],
    [0.1, 0, 0.3, 0.6],
    [0, 0.1, 0.1, 0.8]
])
 
k = 30
d_pi3 = method_multiply(P_pi, 30)
print(f"Multiply {k} times")
print(d_pi3)
 
# obtain the normalized inverse of (I - P_pi)
I = np.eye(P_pi.shape[0])
 
print("Inverse norm direct")
print(inverse_norm_direct(I - P_pi))
 
print("Inverse norm iterative")
print(inverse_norm_iterative(I - P_pi))
 
print("-" * 50)
# discount factor
gamma = 0.9
 
print("Inverse norm direct gamma")
print(inverse_norm_direct(I - gamma * P_pi))
 
print("Inverse norm iterative gamma")
print(inverse_norm_iterative(I - gamma * P_pi))

Multiply 30 times
[[0.03448276 0.10837438 0.13300493 0.72413793] 
 [0.03448276 0.10837438 0.13300493 0.72413793] 
 [0.03448276 0.10837438 0.13300493 0.72413793] 
 [0.03448276 0.10837438 0.13300493 0.72413793]]
Inverse norm direct
[-2.48770265e+16 -2.48770265e+16 -2.48770265e+16 -2.48770265e+16]
[[0.03448276 0.10837438 0.13300493 0.72413793]
 [0.03448276 0.10837438 0.13300493 0.72413793]
 [0.03448276 0.10837438 0.13300493 0.72413793]
 [0.03448276 0.10837438 0.13300493 0.72413793]]
Inverse norm iterative
converged after 100003 iterations
[100002.00000024 100002.00000024 100002.00000024 100002.00000024]
[[0.03449732 0.10837206 0.13301266 0.72411796]
 [0.03448353 0.10838586 0.13300231 0.7241283 ]
 [0.03448353 0.10837157 0.1330166  0.7241283 ]
 [0.03448181 0.10837329 0.13300281 0.72414209]]
--------------------------------------------------
Inverse norm direct gamma
[10. 10. 10. 10.]
[[0.17184995 0.08842402 0.19435892 0.5453671 ]
 [0.04039756 0.21987641 0.10779272 0.63193331]
 [0.04039756 0.08289011 0.24477902 0.63193331]
 [0.02596986 0.09731781 0.11332663 0.7633857 ]]
Inverse norm iterative gamma
converged after 92 iterations
[9.9993144 9.9993144 9.9993144 9.9993144]
[[0.17185937 0.08842265 0.19436313 0.54535485]
 [0.04039797 0.21988405 0.10779099 0.63192699]
 [0.04039797 0.08288836 0.24478668 0.63192699]
 [0.02596928 0.09731705 0.11332528 0.76338839]]
>>> (1 - 0.9**91)/(1-0.9)
9.99931440386759