# 参考 https://qiita.com/python_walker/items/878b72a5215ee424c336
import numpy as np

mat = np.array([
                [1., 5., 1., -1.], 
                [3., 2., 7., 1.], 
                [4., 1., 3., -2.], 
                [1., 6., 4., 3.]
                ])
b = np.array([5., 8., 7., 13.])

for i in range(len(mat)-1):
    for j in range(i+1, len(mat)):
        coef = mat[j][i] / mat[i][i]
        mat[j] -= mat[i] * coef
        b[j] -= b[i] * coef

for i in range(len(mat)-1, 0, -1):
    b[i] /= mat[i][i]
    mat[i] /= mat[i][i]
    for j in range(i):
        b[j] -= b[i] * mat[j][i]
        mat[j][i] = 0

b

array([ 3.5,  1. , -1. ,  2.5])

# 参考 http://d.hatena.ne.jp/aldente39/20110806/1312608456

matrix = np.array([
    [3.0, 2, 7, 1,8],
    [1, 5, 1, -1,5],
    [4, 1, 3, -2,7],
    [1, 6, 4, 3,13]
])

dim=matrix.shape[0]
for i in range(dim):
    num = matrix[i][i]
    for j in range(dim + 1):
        matrix[i][j] = matrix[i][j] / num
    for j in range(dim):
        if i == j:
            pass
        else:
            a = matrix[j][i]
            for k in range(i, dim + 1):
                matrix[j][k] = matrix[j][k]- a * matrix[i][k]
matrix[:,dim]

array([ 3.5,  1. , -1. ,  2.5])

# 参考 https://qiita.com/Dason08/items/4ee5ab1684421f0b6ec9

def solve_CG(A,b,k):
    alpha = 0
    x = np.array([[0],[0],[0],[0]])
    m = A.T * (A*x-b)
    t = -(np.tensordot(m,A.T*(A*x-b)))/(np.tensordot(m,A.T*A*m))
    x = x + t*m

    for i in range(k):
        alpha = - (np.tensordot(m,A.T * A * A.T*(A*x-b)))/(np.tensordot(m, A.T*A*m))
        m = A.T * (A*x-b) + alpha*m
        t = -(np.tensordot(m,A.T*(A*x-b)))/(np.tensordot(m,A.T*A*m))
        x = x + t*m

    return x

A = np.matrix([[3,2,7,1],[1,5,1,-1],[4,1,3,-2],[1,6,4,3]])
b = np.array([[8],[5],[7],[13]])
solve_CG(A,b,3) 

matrix([[ 3.5],
        [ 1. ],
        [-1. ],
        [ 2.5]])