Show how cho_solve gives wrong error and the proper way to do it.

cho_solve_error.py

from scipy.linalg import cho_solve
import numpy as np

def choSolve(L,y,lower=True):
    x = np.copy(y)
    if lower:
        i = 0
        while i < L.shape[0]:
            x[i] /= L[i,i]
            x[i+1:] -= L[i+1:,i]*x[i]
            i += 1       
    return x
y = np.random.uniform(size=3)
a = np.random.uniform(size=[3,3])
a = a.T.dot(a)
L = np.linalg.cholesky(a)
print("Will solve L.x = y")
print("{} . x = {}".format(L,y))
x1 = cho_solve((L,True),y)
print("scipy.linalg.cho_solve result: {}".format(x1))
print("Verify L.x = y: {}".format(np.isclose(L.dot(x1),y)))
x2 = choSolve(L,y)
print("proper results: {}".format(x2))
print("Verify L.x = y: {}".format(np.isclose(L.dot(x2),y)))