Ir al contenido principal

Numerical methods challenge: Day 29

During October (2017) I will write a program per day for some well-known numerical methods in both Python and Julia. It is intended to be an exercise then don't expect the code to be good enough for real use. Also, I should mention that I have almost no experience with Julia, so it probably won't be idiomatic Julia but more Python-like Julia.

Cholesky decomposition

Today we have Cholesky decomposition. It is a factorization of a Hermitian, positive-definite matrix into a lower and upper matrix, the main difference with LU decomposition is that it the lower matrix is the Hermitian transpose of the upper one.

Following are the codes

Python

from __future__ import division, print_function
import numpy as np

def cholesky(mat):
m, _ = mat.shape
mat = mat.copy()
for col in range(m):
print(mat[col, col] -
mat[col, 0:col].dot(mat[col, 0:col]))
mat[col, col] = np.sqrt(mat[col, col] -
mat[col, 0:col].dot(mat[col, 0:col]))
for row in range(col + 1, m):
mat[row, col] = (mat[row, col] -
mat[row, 0:col].dot(mat[col, 0:col]))/mat[col, col]
for row in range(1, m):
mat[0:row, row] = 0
return mat

A = np.array([
[4, -1, 1],
[-1, 4.25, 2.75],
[1, 2.75, 3.5]], dtype=float)
B = cholesky(A)


Julia

function cholesky(mat)
m, _ = size(mat)
mat = copy(mat)
for col = 1:m
mat[col, col] = sqrt(mat[col, col] -
dot(mat[col, 1:col-1], mat[col, 1:col-1]))
for row = col + 1:m
mat[row, col] = (mat[row, col] -
dot(mat[row, 1:col-1], mat[col, 1:col-1]))/mat[col, col]
end
end
for row = 2:m
mat[1:row-1, row] = 0
end
return mat
end

A = [4 -1 1;
-1 4.25 2.75;
1 2.75 3.5]
B = cholesky(A)


As an example we have the matrix

\begin{equation*} A = \begin{bmatrix} 4 &-1 &1\\ -1 &4.25 &2.75\\ 1 &2.75 &3.5 \end{bmatrix} \end{equation*}

And, the answer of both codes is

2.0  0.0  0.0
-0.5  2.0  0.0
0.5  1.5  1.0


Comparison Python/Julia

Regarding number of lines we have: 23 in Python and 22 in Julia. The comparison in execution time is done with %timeit magic command in IPython and @benchmark in Julia.

For Python:

%timeit cholesky(np.eye(100))


with result

100 loops, best of 3: 13 ms per loop


For Julia:

@benchmark cholesky(eye(100))


with result

BenchmarkTools.Trial:
memory estimate:  4.01 MiB
allocs estimate:  20303
--------------
minimum time:     1.010 ms (0.00% GC)
median time:      1.136 ms (0.00% GC)
mean time:        1.370 ms (17.85% GC)
maximum time:     4.652 ms (40.37% GC)
--------------
samples:          3637
evals/sample:     1


In this case, we can say that the Python code is roughly 10 times slower than Julia code.

Comentarios

Comments powered by Disqus