ヒルベルト行列
ヒルベルト行列とその逆行列を求め，両者の行列乗算で単位行列になることを，4 種の言語でやってみた。
Octave と Julia はさすが。
R はちょっと誤差が大きい?
Python は R より誤差が大きい?
Octave
hilb
Return the Hilbert matrix of order N.
invhilb
Return the inverse of the Hilbert matrix of order N.
hilb(5)
ans =
1.0000 0.5000 0.3333 0.2500 0.2000
0.5000 0.3333 0.2500 0.2000 0.1667
0.3333 0.2500 0.2000 0.1667 0.1429
0.2500 0.2000 0.1667 0.1429 0.1250
0.2000 0.1667 0.1429 0.1250 0.1111
invhilb(5)
ans =
25 -300 1050 -1400 630
-300 4800 -18900 26880 -12600
1050 -18900 79380 -117600 56700
-1400 26880 -117600 179200 -88200
630 -12600 56700 -88200 44100
hilb(5) * invhilb(5)
ans =
1.0000 0 0 0 0
0 1.0000 0 0 0
0 0 1.0000 -0.0000 0
0 0 0 1.0000 0
0 0 0 0 1.0000
Julia
# import Pkg; Pkg.add("Nemo")
using Nemo
M = MatrixSpace(QQ, 5, 5)
h = hilbert(M)
Welcome to Nemo version 0.32.2
Nemo comes with absolutely no warranty whatsoever
[ 1 1//2 1//3 1//4 1//5]
[1//2 1//3 1//4 1//5 1//6]
[1//3 1//4 1//5 1//6 1//7]
[1//4 1//5 1//6 1//7 1//8]
[1//5 1//6 1//7 1//8 1//9]
inv(h)
[ 25 -300 1050 -1400 630]
[ -300 4800 -18900 26880 -12600]
[ 1050 -18900 79380 -117600 56700]
[-1400 26880 -117600 179200 -88200]
[ 630 -12600 56700 -88200 44100]
h * inv(h)
[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 1 0]
[0 0 0 0 1]
R
# install.packages("matrixcalc")
library(matrixcalc)
h = hilbert.matrix(5)
h
A matrix: 5 × 5 of type dbl
1.0000000 0.5000000 0.3333333 0.2500000 0.2000000
0.5000000 0.3333333 0.2500000 0.2000000 0.1666667
0.3333333 0.2500000 0.2000000 0.1666667 0.1428571
0.2500000 0.2000000 0.1666667 0.1428571 0.1250000
0.2000000 0.1666667 0.1428571 0.1250000 0.1111111
solve(h)
A matrix: 5 × 5 of type dbl
25 -300 1050 -1400 630
-300 4800 -18900 26880 -12600
1050 -18900 79380 -117600 56700
-1400 26880 -117600 179200 -88200
630 -12600 56700 -88200 44100
h %*% solve(h)
A matrix: 5 × 5 of type dbl
1.000000e+00 -4.547474e-13 -1.818989e-12 0.000000e+00 0.000000e+00
-1.421085e-14 1.000000e+00 0.000000e+00 -5.456968e-12 -2.728484e-12
0.000000e+00 -4.547474e-13 1.000000e+00 0.000000e+00 -9.094947e-13
0.000000e+00 2.273737e-13 0.000000e+00 1.000000e+00 -9.094947e-13
1.421085e-14 -2.273737e-13 0.000000e+00 -1.818989e-12 1.000000e+00
round(h %*% solve(h), 10)
A matrix: 5 × 5 of type dbl
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
Python
from scipy.linalg import hilbert, invhilbert
hilbert(5)
array([[1. , 0.5 , 0.33333333, 0.25 , 0.2 ],
[0.5 , 0.33333333, 0.25 , 0.2 , 0.16666667],
[0.33333333, 0.25 , 0.2 , 0.16666667, 0.14285714],
[0.25 , 0.2 , 0.16666667, 0.14285714, 0.125 ],
[0.2 , 0.16666667, 0.14285714, 0.125 , 0.11111111]])
invhilbert(5)
array([[ 2.500e+01, -3.000e+02, 1.050e+03, -1.400e+03, 6.300e+02],
[-3.000e+02, 4.800e+03, -1.890e+04, 2.688e+04, -1.260e+04],
[ 1.050e+03, -1.890e+04, 7.938e+04, -1.176e+05, 5.670e+04],
[-1.400e+03, 2.688e+04, -1.176e+05, 1.792e+05, -8.820e+04],
[ 6.300e+02, -1.260e+04, 5.670e+04, -8.820e+04, 4.410e+04]])
hilbert(5) @ invhilbert(5)
array([[ 1.00000000e+00, -1.39888101e-13, 6.29496455e-13, 2.65876210e-12, -1.32938105e-12],
[-2.00395256e-14, 1.00000000e+00, -2.34356978e-12, 2.63500333e-12, -1.31750166e-12],
[ 2.34257058e-14, -1.27453603e-13, 1.00000000e+00, -2.93853830e-12, 5.59774449e-13],
[ 0.00000000e+00, -2.27373675e-13, 9.09494702e-13, 1.00000000e+00, 0.00000000e+00],
[-1.80966353e-14, 3.05089287e-13, -3.49720253e-13, 2.36299869e-12, 1.00000000e+00]])
import numpy as np
np.round(hilbert(5) @ invhilbert(5), 10)
array([[ 1., -0., 0., 0., -0.],
[-0., 1., -0., 0., -0.],
[ 0., -0., 1., -0., 0.],
[ 0., -0., 0., 1., 0.],
[-0., 0., -0., 0., 1.]])
