Permalink
Browse files

ENH: linalg: Add a function to compute a Hilbert matrix.

  • Loading branch information...
1 parent 0066902 commit bf05e8d3412e2d96ab6735b6935e2a65f2620c30 Warren Weckesser committed with pv Jun 8, 2011
@@ -53,6 +53,12 @@ A sort keyword has been added to the Schur decomposition routine
(``scipy.linalg.schur``) to allow the sorting of eigenvalues in
the resultant Schur form.
+Additional special matrices (``scipy.linalg``)
+----------------------------------------------
+
+The functions ``hilbert`` and ``invhilbert`` were added to ``scipy.linalg``.
+
+
Deprecated features
===================
@@ -88,6 +88,8 @@ Special Matrices
companion
hadamard
hankel
+ hilbert
+ invhilbert
kron
leslie
toeplitz
@@ -843,6 +843,10 @@ Hadamard `scipy.linalg.hadamard` Construct a Hadamard matrix.
-------------------- ------------------------- ---------------------------------------------------------
Hankel `scipy.linalg.hankel` Construct a Hankel matrix.
-------------------- ------------------------- ---------------------------------------------------------
+Hilbert `scipy.linalg.hilbert` Construct a Hilbert matrix.
+-------------------- ------------------------- ---------------------------------------------------------
+Inverse Hilbert `scipy.linalg.invhilbert` Construct the inverse of a Hilbert matrix.
+-------------------- ------------------------- ---------------------------------------------------------
Leslie `scipy.linalg.leslie` Create a Leslie matrix.
-------------------- ------------------------- ---------------------------------------------------------
Toeplitz `scipy.linalg.toeplitz` Construct a Toeplitz matrix.
@@ -85,6 +85,8 @@
companion - Companion matrix
hadamard - Hadamard matrix of order 2**n
hankel - Hankel matrix
+ hilbert - Hilbert matrix
+ invhilbert - Inverse Hilbert matrix
kron - Kronecker product of two arrays
leslie - Leslie matrix
toeplitz - Toeplitz matrix
@@ -1,9 +1,11 @@
import math
import numpy as np
+from scipy.misc import comb
__all__ = ['tri', 'tril', 'triu', 'toeplitz', 'circulant', 'hankel',
- 'hadamard', 'leslie', 'all_mat', 'kron', 'block_diag', 'companion']
+ 'hadamard', 'leslie', 'all_mat', 'kron', 'block_diag', 'companion',
+ 'hilbert', 'invhilbert']
#-----------------------------------------------------------------------------
@@ -571,3 +573,99 @@ def companion(a):
c[0] = first_row
c[range(1,n-1), range(0, n-2)] = 1
return c
+
+
+def hilbert(n):
+ """Create a Hilbert matrix of order n.
+
+ Returns the `n` by `n` array with entries `h[i,j] = 1 / (i + j + 1)`.
+
+ Parameters
+ ----------
+ n : int
+ The size of the array to create.
+
+ Returns
+ -------
+ h : ndarray with shape (n, n)
+ The Hilber matrix.
+
+ Notes
+ -----
+ .. versionadded:: 0.10.0
+
+ Examples
+ --------
+ >>> hilbert(3)
+ array([[ 1. , 0.5 , 0.33333333],
+ [ 0.5 , 0.33333333, 0.25 ],
+ [ 0.33333333, 0.25 , 0.2 ]])
+
+ """
+ values = 1.0 / (1.0 + np.arange(2 * n - 1))
+ h = hankel(values[:n], r=values[n-1:])
+ return h
+
+
+def invhilbert(n, exact=False):
+ """Compute the inverse of the Hilbert matrix of order `n`.
+
+ Parameters
+ ----------
+ n : int
+ The order of the Hilbert matrix.
+ exact : bool
+ If False, the data type of the array that is returned in np.float64,
+ and the array is an approximation of the inverse.
+ If True, the array is exact integer array. To represent the exact
+ inverse when n > 14, the returned array is an object array of long
+ integers. For n <= 14, the exact inverse is returned as an array
+ with data type np.int64.
+
+ Returns
+ -------
+ invh : ndarray with shape (n, n)
+ The data type of the array is np.float64 is exact is False.
+ If exact is True, the data type is either np.int64 (for n <= 14)
+ or object (for n > 14). In the latter case, the objects in the
+ array will be long integers.
+
+ Notes
+ -----
+ .. versionadded:: 0.10.0
+
+ Examples
+ --------
+ >>> invhilbert(4)
+ array([[ 16., -120., 240., -140.],
+ [ -120., 1200., -2700., 1680.],
+ [ 240., -2700., 6480., -4200.],
+ [ -140., 1680., -4200., 2800.]])
+ >>> invhilbert(4, exact=True)
+ array([[ 16, -120, 240, -140],
+ [ -120, 1200, -2700, 1680],
+ [ 240, -2700, 6480, -4200],
+ [ -140, 1680, -4200, 2800]], dtype=int64)
+ >>> invhilbert(16)[7,7]
+ 4.2475099528537506e+19
+ >>> invhilbert(16, exact=True)[7,7]
+ 42475099528537378560L
+ """
+ if exact:
+ if n > 14:
+ dtype = object
+ else:
+ dtype = np.int64
+ else:
+ dtype = np.float64
+ invh = np.empty((n, n), dtype=dtype)
+ for i in xrange(n):
+ for j in xrange(0, i + 1):
+ s = i + j
+ invh[i, j] = ((-1)**s * (s + 1) *
+ comb(n + i, n - j - 1, exact) *
+ comb(n + j, n - i - 1, exact) *
+ comb(s, i, exact) ** 2)
+ if i != j:
+ invh[j, i] = invh[i, j]
+ return invh
@@ -2,10 +2,13 @@
from numpy import arange, add, array, eye, copy
from numpy.testing import TestCase, run_module_suite, assert_raises, \
- assert_equal, assert_array_equal
+ assert_equal, assert_array_equal, assert_array_almost_equal, \
+ assert_allclose
from scipy.linalg import toeplitz, hankel, circulant, hadamard, leslie, \
- companion, tri, triu, tril, kron, block_diag
+ companion, tri, triu, tril, kron, block_diag, \
+ hilbert, invhilbert
+from numpy.linalg import cond
def get_mat(n):
@@ -265,5 +268,174 @@ def test_basic(self):
assert_array_equal(a, expected)
+class TestHilbert(TestCase):
+
+ def test_basic(self):
+ h3 = array([[1.0, 1/2., 1/3.],
+ [1/2., 1/3., 1/4.],
+ [1/3., 1/4., 1/5.]])
+ assert_array_almost_equal(hilbert(3), h3)
+
+ assert_array_equal(hilbert(1), [[1.0]])
+
+ h0 = hilbert(0)
+ assert_equal(h0.shape, (0,0))
+
+
+class TestInvHilbert(TestCase):
+
+ def test_basic(self):
+ invh1 = array([[1]])
+ assert_array_equal(invhilbert(1, exact=True), invh1)
+ assert_array_equal(invhilbert(1), invh1)
+
+ invh2 = array([[ 4, -6],
+ [-6, 12]])
+ assert_array_equal(invhilbert(2, exact=True), invh2)
+ assert_array_almost_equal(invhilbert(2), invh2)
+
+ invh3 = array([[ 9, -36, 30],
+ [-36, 192, -180],
+ [30, -180, 180]])
+ assert_array_equal(invhilbert(3, exact=True), invh3)
+ assert_array_almost_equal(invhilbert(3), invh3)
+
+ invh4 = array([[ 16, -120, 240, -140],
+ [-120, 1200, -2700, 1680],
+ [ 240, -2700, 6480, -4200],
+ [-140, 1680, -4200, 2800]])
+ assert_array_equal(invhilbert(4, exact=True), invh4)
+ assert_array_almost_equal(invhilbert(4), invh4)
+
+ invh5 = array([[ 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]])
+ assert_array_equal(invhilbert(5, exact=True), invh5)
+ assert_array_almost_equal(invhilbert(5), invh5)
+
+ invh17 = array([
+ [289, -41616, 1976760, -46124400, 629598060, -5540462928,
+ 33374693352, -143034400080, 446982500250, -1033026222800,
+ 1774926873720, -2258997839280, 2099709530100, -1384423866000,
+ 613101997800, -163493866080, 19835652870],
+ [-41616, 7990272, -426980160, 10627061760, -151103534400, 1367702848512,
+ -8410422724704, 36616806420480, -115857864064800, 270465047424000,
+ -468580694662080, 600545887119360, -561522320049600, 372133135180800,
+ -165537539406000, 44316454993920, -5395297580640],
+ [1976760, -426980160, 24337869120, -630981792000, 9228108708000,
+ -85267724461920, 532660105897920, -2348052711713280, 7504429831470000,
+ -17664748409880000, 30818191841236800, -39732544853164800,
+ 37341234283298400, -24857330514030000, 11100752642520000,
+ -2982128117299200, 364182586693200],
+ [-46124400, 10627061760, -630981792000, 16826181120000,
+ -251209625940000, 2358021022156800, -14914482965141760,
+ 66409571644416000, -214015221119700000, 507295338950400000,
+ -890303319857952000, 1153715376477081600, -1089119333262870000,
+ 727848632044800000, -326170262829600000, 87894302404608000,
+ -10763618673376800],
+ [629598060, -151103534400, 9228108708000,
+ -251209625940000, 3810012660090000, -36210360321495360,
+ 231343968720664800, -1038687206500944000, 3370739732635275000,
+ -8037460526495400000, 14178080368737885600, -18454939322943942000,
+ 17489975175339030000, -11728977435138600000, 5272370630081100000,
+ -1424711708039692800, 174908803442373000],
+ [-5540462928, 1367702848512, -85267724461920, 2358021022156800,
+ -36210360321495360, 347619459086355456, -2239409617216035264,
+ 10124803292907663360, -33052510749726468000, 79217210949138662400,
+ -140362995650505067440, 183420385176741672960, -174433352415381259200,
+ 117339159519533952000, -52892422160973595200, 14328529177999196160,
+ -1763080738699119840],
+ [33374693352, -8410422724704, 532660105897920,
+ -14914482965141760, 231343968720664800, -2239409617216035264,
+ 14527452132196331328, -66072377044391477760, 216799987176909536400,
+ -521925895055522958000, 928414062734059661760, -1217424500995626443520,
+ 1161358898976091015200, -783401860847777371200, 354015418167362952000,
+ -96120549902411274240, 11851820521255194480],
+ [-143034400080, 36616806420480, -2348052711713280, 66409571644416000,
+ -1038687206500944000, 10124803292907663360, -66072377044391477760,
+ 302045152202932469760, -995510145200094810000, 2405996923185123840000,
+ -4294704507885446054400, 5649058909023744614400,
+ -5403874060541811254400, 3654352703663101440000,
+ -1655137020003255360000, 450325202737117593600, -55630994283442749600],
+ [446982500250, -115857864064800, 7504429831470000, -214015221119700000,
+ 3370739732635275000, -33052510749726468000, 216799987176909536400,
+ -995510145200094810000, 3293967392206196062500,
+ -7988661659013106500000, 14303908928401362270000,
+ -18866974090684772052000, 18093328327706957325000,
+ -12263364009096700500000, 5565847995255512250000,
+ -1517208935002984080000, 187754605706619279900],
+ [-1033026222800, 270465047424000, -17664748409880000,
+ 507295338950400000, -8037460526495400000, 79217210949138662400,
+ -521925895055522958000, 2405996923185123840000,
+ -7988661659013106500000, 19434404971634224000000,
+ -34894474126569249192000, 46141453390504792320000,
+ -44349976506971935800000, 30121928988527376000000,
+ -13697025107665828500000, 3740200989399948902400,
+ -463591619028689580000],
+ [1774926873720, -468580694662080,
+ 30818191841236800, -890303319857952000, 14178080368737885600,
+ -140362995650505067440, 928414062734059661760, -4294704507885446054400,
+ 14303908928401362270000, -34894474126569249192000,
+ 62810053427824648545600, -83243376594051600326400,
+ 80177044485212743068000, -54558343880470209780000,
+ 24851882355348879230400, -6797096028813368678400, 843736746632215035600],
+ [-2258997839280, 600545887119360, -39732544853164800,
+ 1153715376477081600, -18454939322943942000, 183420385176741672960,
+ -1217424500995626443520, 5649058909023744614400,
+ -18866974090684772052000, 46141453390504792320000,
+ -83243376594051600326400, 110552468520163390156800,
+ -106681852579497947388000, 72720410752415168870400,
+ -33177973900974346080000, 9087761081682520473600,
+ -1129631016152221783200],
+ [2099709530100, -561522320049600, 37341234283298400,
+ -1089119333262870000, 17489975175339030000, -174433352415381259200,
+ 1161358898976091015200, -5403874060541811254400,
+ 18093328327706957325000, -44349976506971935800000,
+ 80177044485212743068000, -106681852579497947388000,
+ 103125790826848015808400, -70409051543137015800000,
+ 32171029219823375700000, -8824053728865840192000,
+ 1098252376814660067000],
+ [-1384423866000, 372133135180800,
+ -24857330514030000, 727848632044800000, -11728977435138600000,
+ 117339159519533952000, -783401860847777371200, 3654352703663101440000,
+ -12263364009096700500000, 30121928988527376000000,
+ -54558343880470209780000, 72720410752415168870400,
+ -70409051543137015800000, 48142941226076592000000,
+ -22027500987368499000000, 6049545098753157120000,
+ -753830033789944188000],
+ [613101997800, -165537539406000,
+ 11100752642520000, -326170262829600000, 5272370630081100000,
+ -52892422160973595200, 354015418167362952000, -1655137020003255360000,
+ 5565847995255512250000, -13697025107665828500000,
+ 24851882355348879230400, -33177973900974346080000,
+ 32171029219823375700000, -22027500987368499000000,
+ 10091416708498869000000, -2774765838662800128000, 346146444087219270000],
+ [-163493866080, 44316454993920, -2982128117299200, 87894302404608000,
+ -1424711708039692800, 14328529177999196160, -96120549902411274240,
+ 450325202737117593600, -1517208935002984080000, 3740200989399948902400,
+ -6797096028813368678400, 9087761081682520473600,
+ -8824053728865840192000, 6049545098753157120000,
+ -2774765838662800128000, 763806510427609497600, -95382575704033754400],
+ [19835652870, -5395297580640, 364182586693200, -10763618673376800,
+ 174908803442373000, -1763080738699119840, 11851820521255194480,
+ -55630994283442749600, 187754605706619279900, -463591619028689580000,
+ 843736746632215035600, -1129631016152221783200, 1098252376814660067000,
+ -753830033789944188000, 346146444087219270000, -95382575704033754400,
+ 11922821963004219300]
+ ])
+ assert_array_equal(invhilbert(17, exact=True), invh17)
+ assert_allclose(invhilbert(17), invh17.astype(float), rtol=1e-12)
+
+ def test_inverse(self):
+ for n in xrange(1, 10):
+ a = hilbert(n)
+ b = invhilbert(n)
+ # The Hilbert matrix is increasingly badly conditioned,
+ # so take that into account in the test
+ c = cond(a)
+ assert_allclose(a.dot(b), eye(n), atol=1e-15*c, rtol=1e-15*c)
+
if __name__ == "__main__":
run_module_suite()

0 comments on commit bf05e8d

Please sign in to comment.