Merge pull request #228 from glorotxa/SpectralRadiusBound_review

theano/sandbox/linalg/ops.py

@@ -691,7 +691,14 @@ def spectral_radius_bound(X, log2_exponent):
 
     From V.Pan, 1990. "Estimating the Extremal Eigenvalues of a Symmetric
     Matrix", Computers Math Applic. Vol 20 n. 2 pp 17-22.
+    Rq: an efficient algorithm, not used here, is defined in this paper.
     """
+    if X.type.ndim != 2:
+        raise TypeError('spectral_radius_bound requires a matrix argument', X)
+    if not(isintance(log2_exponent,int)):
+        raise TypeError('spectral_radius_bound requires a integer exponent', log2_exponent)
+    if log2_exponent<=0:
+        raise ValueError('spectral_radius_bound requires a strictly positive exponent', log2_exponent)
     XX = X
     for i in xrange(log2_exponent):
         XX = tensor.dot(XX, XX)

theano/sandbox/linalg/tests/test_linalg.py

@@ -23,7 +23,7 @@ from theano.sandbox.linalg.ops import (cholesky,
                                        #PSD_hint,
                                        trace,
                                        matrix_dot,
-                                       #spectral_radius_bound
+                                       spectral_radius_bound
                                        )
 
 from nose.plugins.skip import SkipTest
@@ -281,4 +281,48 @@ def test_trace():
         ok = True
     assert ok
 
-
+def test_spectral_radius_bound():
+    tol = 10**(-6)
+    rng = numpy.random.RandomState(utt.fetch_seed())
+    x = theano.tensor.matrix()
+    radius_bound = spectral_radius_bound(x, 5)
+    f = theano.function([x], radius_bound)
+    
+    shp = (3, 3)
+    m = rng.rand(*shp).astype(config.floatX)
+    m = numpy.cov(m)
+    radius_bound_theano = f(m)
+    
+    # test the approximation
+    mm = m
+    for i in range(5):
+        mm = numpy.dot(mm, mm)
+    radius_bound_numpy = numpy.trace(mm)**(2**(-5))
+    assert abs(radius_bound_numpy - radius_bound_theano) < tol
+    
+    # test the bound 
+    eigen_val = numpy.linalg.eig(m)
+    assert (eigen_val[0].max() - radius_bound_theano) < tol
+    
+    # test type errors
+    xx = theano.tensor.vector()
+    ok = False
+    try:
+        spectral_radius_bound(xx, 5)
+    except TypeError:
+        ok = True
+    assert ok
+    ok = False
+    try:
+        spectral_radius_bound(x, 5.)
+    except TypeError:
+        ok = True
+    assert ok
+    
+    # test value error
+    ok = False
+    try:
+        spectral_radius_bound(x, -5)
+    except ValueError:
+        ok = True
+    assert ok