C++ vs Python vs JavaScript vs Julia

今回は、Gaussの消去法を使って、C++、Python、JavaScript、Juliaの計算時間を計測しました。
Pythonに関しては、JITコンパイラであるNumba無しと有りでどうなるかを評価しました。
また、Numbaで@numba.jitを付けた関数を２回以上呼び出した場合、１回目と２回目以降では計算速度が異なります。２回目以降はcacheから関数を呼び出すため高速になるので、これについても評価しました。

〈テスト環境〉
・CPU：Intel Celeron Processor 1005M（1.9GHz）
・OS：Lubuntu 20.04.4
・Python 3.8.10
・Numpy 1.21.1
・Numba 0.55.1
・llvmlite 0.38.0

結果は表１のようになりました。
PythonはNumba無しからNumba有りにすると、１回目の計算では350倍、２回目の計算では570倍高速化しました。
また、JavaScript（V8）が意外と健闘しました。Windows 8.1ではJuliaよりもわずかに速かったです。V8もNumbaと同様に２回目以降の計算を高速化してくれると嬉しいのですが。。。
Juliaに関してはC++と同レベルだと思っていただけにちょっと残念な結果に終わりましたが、ベクトルや行列の演算ができるので、今後AIの分野での利用が拡大するのではないかと思います。

注）Pythonでクラスの中の関数（メソッド）に@numba.jitを付けるとWarningが表示され高速化できなかったです。
注）Juliaの計算では２次元配列 Array{Float64, 2}:a[i, j] ではなく、配列の配列 Vector{Vector{Float64}}:a[i][j] を用いています。２次元配列を用いた場合、計算時間が6.888756secとなり、8倍以上遅くなりました。
？）gauss.pyのプログラム（Numba有り）で、set(a, b)、solve(a, b) を set(n, a, b)、solve(n, a, b) にすると十数％遅くなりました。
また、#n = len(a) で # を外しても十数％遅くなりました。どうしてこんなに遅くなるのでしょう？

表１．計算時間（Gaussの消去法）

プログラミング言語	計算時間	C++を1としたときの計算時間
C++（Clang）	0.508967 sec	1.000
Python（Numba無し１回目の計算）	377.815687 sec	742.319
Python（Numba無し２回目の計算）	378.716451 sec	744.088
Python（Numba有り１回目の計算）	1.074311 sec	2.111
Python（Numba有り２回目の計算）	0.655199 sec	1.287
JavaScript（Node.js:V8）	0.998315 sec	1.961
Julia	0.796610 sec	1.565

Gaussの消去法は連立一次方程式の解法の一つです。今回のプログラムではピボットを選択しない方法を採用しているので、対角要素に０がある場合は使えません。 プログラムでは、係数行列は対角要素より下の要素は０、それ以外の要素は１、非同次項は全ての要素を１に設定しています。
     a x = b
     a：係数行列（n×nの正方行列）
     b：非同次項（n次元の列ベクトル）
     x：求める解（n次元の列ベクトル）

今回、計算時間を計測するに当たって用いたプログラムは以下のようになります。

gauss.cpp（Clang）

/*
$ sudo apt update
$ sudo apt install clang -y
$ clang++ gauss.cpp -o gauss -O3
$ ./gauss
*/

#include <chrono>
#include <cstdio>

class Gauss {
  static const int n = 1000;
  double a[n][n] = {{}};
  double b[n] = {};
public:
  auto set(double a, double b) -> void {
    for (int i = 0; i < n; i++) {
      for (int j = i; j < n; j++) {
        this->a[i][j] = a;
      }
      this->b[i] = b;
    }
  }
  auto solve() -> void {
    for (int k = 0; k < n - 1; k++) {
      for (int i = k + 1; i < n; i++) {
        const double s = a[i][k] / a[k][k];
        for (int j = k + 1; j < n; j++) {
          a[i][j] -= s * a[k][j];
        }
        b[i] -= s * b[k];
      }
    }
    for (int i = n - 1; i >= 0; i--) {
      for (int j = i + 1; j < n; j++) {
        b[i] -= a[i][j] * b[j];
      }
      b[i] /= a[i][i];
    }
  }
  auto print() -> void {
    for (int i = 0; i < n; i++) {
      printf("x[%d] = %.10E\n", i, b[i]);
    }
  }
};

auto main() -> int {
  const auto gauss = new Gauss;
  gauss->set(1, 1);
  const auto start = std::chrono::system_clock::now();
  gauss->solve();
  const auto stop = std::chrono::system_clock::now();
  gauss->print();
  printf("%.6f%s\n", 0.000001 * std::chrono::duration_cast<std::chrono::microseconds>(stop - start).count(), "sec");
  delete gauss;
}
 

下記のNumba有り（JIT有り）のプログラムで@numba.jitをコメントアウトすれば、Numba無し（JIT無し）のプログラムになります。
また、@numba.jit(cache=True)の方を有効にすると、２回目のプログラムの実行から１回目以降の計算をcacheから呼び出して高速化します。

gauss.py（Python×Numba）

"""
$ sudo apt update
$ sudo apt install python3 python3-pip -y
$ pip3 install numba==0.55.1 numpy==1.21.1
$ python3 gauss.py
"""
import numba
import numpy as np
import time

@numba.jit
#@numba.jit(cache=True)
def set(a, b, a0, b0):
  #n = len(a)
  for i in range(n):
    for j in range(i, n):
      a[i][j] = a0
    b[i] = b0

@numba.jit
#@numba.jit(cache=True)
def solve(a, b):
  #n = len(a)
  for k in range(n - 1):
    for i in range(k + 1, n):
      s = a[i][k] / a[k][k]
      for j in range(k + 1, n):
        a[i][j] -= s * a[k][j]
      b[i] -= s * b[k]
  for i in range(n - 1, - 1, - 1):
    for j in range(i + 1, n):
      b[i] -= a[i][j] * b[j]
    b[i] /= a[i][i]

def printf(b):
  for i in range(len(b)):
    print(f"x[{i}] = {b[i]:.10E}")

#１回目の計算
n = 1000
a = np.zeros((n, n))
b = np.zeros(n)
set(a, b, 1, 1)
start = time.time()
solve(a, b)
stop = time.time()
#printf(b)
print(f"{(stop - start):.6f}sec")

#２回目の計算
n = 1000
a = np.zeros((n, n))
b = np.zeros(n)
set(a, b, 1, 1)
start = time.time()
solve(a, b)
stop = time.time()
#printf(b)
print(f"{(stop - start):.6f}sec")
 

gauss.js（Node.js, Bun）

/*
$ sudo apt update
$ sudo apt install npm -y
$ sudo npm install n -g
$ sudo n stable
$ node gauss.js
*/
/*
$ sudo apt update
$ sudo apt install curl -y
$ curl https://bun.sh/install | bash
$ export PATH=~/.bun/bin:$PATH
$ bun gauss.js
*/

"use strict";
const performance = require("perf_hooks").performance;

class Gauss {
  constructor(n) {
    this.a = Array(n).fill().map(() => new Float64Array(n).fill(0));
    this.b = new Float64Array(n).fill(0);
  }
  set(a, b) {
    const n = this.a.length;
    for (let i = 0; i < n; i++) {
      for (let j = i; j < n; j++) {
        this.a[i][j] = a;
      }
      this.b[i] = b;
    }
  }
  solve() {
    const n = this.a.length;
    for (let k = 0; k < n - 1; k++) {
      for (let i = k + 1; i < n; i++) {
        const s = this.a[i][k] / this.a[k][k];
        for (let j = k + 1; j < n; j++) {
          this.a[i][j] -= s * this.a[k][j];
        }
        this.b[i] -= s * this.b[k];
      }
    }
    for (let i = n - 1; i >= 0; i--) {
      for (let j = i + 1; j < n; j++) {
        this.b[i] -= this.a[i][j] * this.b[j];
      }
      this.b[i] /= this.a[i][i];
    }
  }
  print() {
    for (const i in this.b) {
      console.log(`x[${i}] = ${this.b[i].toExponential(10)}`);
    }
  }
}

const gauss = new Gauss(1000);
gauss.set(1, 1);
const start = performance.now();
gauss.solve();
const stop = performance.now();
gauss.print();
console.log(`${(0.001 * (stop - start)).toFixed(6)}sec`);
 

Juliaにはclassがないのですが、structを使ってオブジェクト指向風に書いてみました。 上記のJavaScriptのプログラムと比較すると、オブジェクトobjに付随するメソッドmの呼び出しが、構造体strを引数に持つ関数fの呼び出しに変わったのがわかると思います。
   obj.m(args)　→　f(str, args)

gauss.jl（Julia）

#=
$ sudo apt update
$ sudo apt install julia -y
$ julia gauss.jl
=#

using Printf

struct Gauss
  a::Vector{Vector{Float64}}
  b::Vector{Float64}
end

function Gauss(n::Int)
  a = [zeros(Float64, n) for i = 1 : n]
  b = zeros(Float64, n)
  Gauss(a, b)
end

function set!(gauss::Gauss, a::Float64, b::Float64)
  n = length(gauss.a)
  for i = 1 : n
    for j = i : n
      gauss.a[i][j] = a
    end
    gauss.b[i] = b
  end
end

function solve!(gauss::Gauss)
  n = length(gauss.a)
  for k = 1 : n - 1
    for i = k + 1 : n
      s = gauss.a[i][k] / gauss.a[k][k]
      for j = k + 1 : n
        gauss.a[i][j] -= s * gauss.a[k][j]
      end
      gauss.b[i] -= s * gauss.b[k]
    end
  end
  for i = n : - 1 : 1
    for j = i + 1 : n
      gauss.b[i] -= gauss.a[i][j] * gauss.b[j]
    end
    gauss.b[i] /= gauss.a[i][i]
  end
end

function print(gauss::Gauss)
  for i = 1 : length(gauss.b)
    @printf("x[%d] = %.10E\n", i, gauss.b[i])
  end
end

const gauss = Gauss(1000)
set!(gauss, 1.0, 1.0)
@time begin
  solve!(gauss)
end
#print(gauss)
 

注）Vector{Float64}はArray{Float64, 1}と同じです。
注）Juliaでは引数に渡された変数を変更する関数（破壊的な関数）には関数名の最後に！を付けるそうですが、付けなくてもエラーにはなりません。

♦処理系のインストールと実行

今回は、Lubuntu 20.04.4 LTS（Focal Fossa）でテストしました。Lubuntuは、ダウンロードサイトからISOファイルをダウンロードし、ISOファイルを開き、その中身を全てフォーマット済みのUSBメモリにコピーして、F12キーを連打しながら再起動（途中でUSBメモリを挿入したドライブを選択）すれば使えるようになります。 この状態でソフトウェアのインストールが可能です。但し、ディスクではなくメモリにインストールされるだけなので、シャットダウンすれば消えてしまいます。

 

〈Clangのインストールと実行〉

$ sudo apt update
$ sudo apt install clang -y
$ clang++ gauss.cpp -o gauss -O3
$ ./gauss

〈PythonとNumbaのインストールと実行〉

$ sudo apt update
$ sudo apt install python3 python3-pip -y
$ pip3 install numba==0.55.1 numpy==1.21.1
$ python3 gauss.py

注）Numpyはバージョンによってインストール時にエラーになることがあります。numpy-1.21.1はOKでした。

 

〈Node.jsのインストールと実行〉

$ sudo apt update
$ sudo apt install npm -y
$ sudo npm install n -g
$ sudo n stable
$ node gauss.js

〈Bunのインストールと実行〉

$ sudo apt update
$ sudo apt install curl -y
$ curl https://bun.sh/install | bash
$ export PATH=~/.bun/bin:$PATH
$ bun gauss.js

〈Juliaのインストールと実行〉

$ sudo apt update
$ sudo apt install julia -y
$ julia gauss.jl

参考記事

• Numba
• NumbaでざっくりPython高速化
• NumbaとCythonの実行速度を比較してみた

• Julia
• Julia docs：Constructors
• Juliaで行列を高速に更新する