add FDA for comparison

Author: rflamary

examples/plot_WDA.py

@@ -11,7 +11,7 @@
 import matplotlib.pylab as pl
 import ot
 from ot.datasets import get_1D_gauss as gauss
-from ot.dr import wda
+from ot.dr import wda, fda
 
 
 #%% parameters
@@ -36,7 +36,12 @@
 pl.title('Discriminant dimensions')
 
 
-#%% plot distributions and loss matrix
+#%% Comlpute FDA
+p=2
+
+Pfda,projfda = fda(xs,ys,p)
+
+#%% Compute WDA
 p=2
 reg=1
 k=10
@@ -46,8 +51,8 @@
 
 #%% plot samples
 
-xsp=proj(xs)
-xtp=proj(xt)
+xsp=projfda(xs)
+xtp=projfda(xt)
 
 pl.figure(1,(10,5))
 

ot/dr.py

@@ -7,6 +7,7 @@
 from pymanopt.manifolds import Stiefel
 from pymanopt import Problem
 from pymanopt.solvers import SteepestDescent, TrustRegions
+import scipy.linalg as la
 
 def dist(x1,x2):
     """ Compute squared euclidean distance between samples (autograd)
@@ -32,9 +33,73 @@ def split_classes(X,y):
     """
     lstsclass=np.unique(y)
     return [X[y==i,:].astype(np.float32) for i in lstsclass]
+
+
+def fda(X,y,p=2,reg=1e-16):
+    """ 
+    Fisher Discriminant Analysis 
+      
     
+    Parameters
+    ----------
+    X : numpy.ndarray (n,d)
+        Training samples
+    y : np.ndarray (n,)
+        labels for training samples
+    p : int, optional
+        size of dimensionnality reduction
+    reg : float, optional
+        Regularization term >0 (ridge regularization)
 
 
+    Returns
+    -------
+    P : (d x p) ndarray
+        Optimal transportation matrix for the given parameters
+    proj : fun
+        projection function including mean centering
+
+    
+    """    
+    
+    mx=np.mean(X)
+    X-=mx.reshape((1,-1))
+    
+    # data split between classes
+    d=X.shape[1]
+    xc=split_classes(X,y)  
+    nc=len(xc)
+    
+    p=min(nc-1,p)
+    
+    Cw=0
+    for x in xc:
+        Cw+=np.cov(x,rowvar=False)
+    Cw/=nc
+    
+    mxc=np.zeros((d,nc))
+    
+    for i in range(nc):
+        mxc[:,i]=np.mean(xc[i])
+        
+    mx0=np.mean(mxc,1)
+    Cb=0
+    for i in range(nc):
+        Cb+=(mxc[:,i]-mx0).reshape((-1,1))*(mxc[:,i]-mx0).reshape((1,-1))
+        
+    w,V=la.eig(Cb,Cw+reg*np.eye(d))
+    
+    idx=np.argsort(w.real)
+    
+    Popt=V[:,idx[-p:]]
+    
+    
+    
+    def proj(X):
+        return (X-mx.reshape((1,-1))).dot(Popt)
+    
+    return Popt, proj
+
 def wda(X,y,p=2,reg=1,k=10,solver = None,maxiter=100,verbose=0):
     """ 
     Wasserstein Discriminant Analysis [11]_
@@ -73,16 +138,13 @@ def wda(X,y,p=2,reg=1,k=10,solver = None,maxiter=100,verbose=0):
     P : (d x p) ndarray
         Optimal transportation matrix for the given parameters
     proj : fun
-        projectiuon function including mean centering
+        projection function including mean centering
 
 
     References
     ----------
 
     .. [11] Flamary, R., Cuturi, M., Courty, N., & Rakotomamonjy, A. (2016). Wasserstein Discriminant Analysis. arXiv preprint arXiv:1608.08063.
-
-    
-    
     
     """
 
@@ -131,3 +193,6 @@ def proj(X):
         return (X-mx.reshape((1,-1))).dot(Popt)
 
     return Popt, proj
+
+
+