Changeset 834

Timestamp:

10/13/15 13:47:42 (7 months ago)

Author:

mmckerns

Message:

extended svm and svc KernelMatrix? to n-dimensions

Location:

mystic/mystic

Files:

• svctools.py (modified) (2 diffs)
• svmtools.py (modified) (3 diffs)

mystic/mystic/svctools.py

@@ -2 +2 @@
 #
 # Author: Patrick Hung (patrickh @caltech)
+# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
 # Copyright (c) 1997-2015 California Institute of Technology.
 # License: 3-clause BSD.  The full license text is available at:
@@ -9 +10 @@
 """
 
-from numpy import zeros, multiply, ndarray, vectorize, array, dot, transpose, diag, sum
+from numpy import multiply, asarray, dot, transpose, sum
 
 def KernelMatrix(X, k=dot):
-    n = X.shape[0]
-    Q = zeros((n,n))
-    for i in range(n):
-       for j in range(i, n):
-           Q[i,j] = k(X[i,:],X[j,:])
-    return Q + transpose(Q) - diag(Q.diagonal())
+    "inner product of X with self, using k as elementwise product function"
+    # the following is tensordot(X,X,axes=(-1,-1)), with dot --> k
+    # 3-clause BSD (see: v1.7.2 http://docs.scipy.org/doc/numpy/license.html)
+    X = asarray(X)
+    Xs = X.shape
+    ndX = len(X.shape)
+    nX = 1
+    axes = [ndX - 1]
+
+    # Move the axes to sum over to the end of "a"
+    # and to the front of "b" in inner(a,b)
+    notin = [_ for _ in range(ndX) if _ not in axes]
+    newaxes_a = notin + axes
+    N2 = 1
+    for axis in axes: N2 *= Xs[axis]
+    newshape_a = (-1, N2)
+    olda = [Xs[axis] for axis in notin]
+
+    newaxes_b = axes + notin
+    N2 = 1
+    for axis in axes: N2 *= Xs[axis]
+    newshape_b = (N2, -1)
+    oldb = [Xs[axis] for axis in notin]
+
+    at = X.transpose(newaxes_a).reshape(newshape_a)
+    bt = X.transpose(newaxes_b).reshape(newshape_b)
+    return k(at,bt).reshape(olda + oldb)
+
 
 def WeightVector(alpha, X, y):

mystic/mystic/svmtools.py

@@ -2 +2 @@
 #
 # Author: Patrick Hung (patrickh @caltech)
+# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
 # Copyright (c) 1997-2015 California Institute of Technology.
 # License: 3-clause BSD.  The full license text is available at:
@@ -9 +10 @@
 """
 
-from numpy import zeros, multiply, ndarray, vectorize, array
+from numpy import multiply, ndarray, vectorize, array, asarray
 
-def InnerProduct(i1, i2):
+def InnerProduct(i1, i2): # numpy.multiply(x,x)
     return i1 * i2
 
@@ -17 +18 @@
     return 1. + i1 * i2
 
-def KernelMatrix(X, k=InnerProduct):
-    n = X.size
-    Q = zeros((n,n))
-    for i in range(n):
-       for j in range(n):
-           # dumb, but how else to do outer products of arbitrary functions
-           # without going through ufunc C-api ?
-           Q[i,j] = k(X[i],X[j])
-    return Q
+def KernelMatrix(X, k=multiply):
+    "outer product of X with self, using k as elementwise product function"
+    X = asarray(X).ravel()
+    return k(X[:,None], X[None,:])
+
 
 def SupportVectors(alpha, eps=0):