opencv2笔记03-线性判别分析（LDA）

1. 线性判别分析（二类情况）

判别准则函数：
$\Large J(w) = \frac { |\tilde{\mu_1} - \tilde{\mu_2}| ^2 }{ \tilde{ S_1} ^2 + \tilde{ S_2} ^2}$
其中：

1.d维x空间
各类样本的均值向量：
$\large \mu_i= \frac{ 1}{N_i } \sum_{x\in \omega_i} x$
类内离散度矩阵：
$\large {S_i}=\sum_{x\in \omega_i} ( x-\mu_i)( x-\mu_i)^T$
总类内离散度矩阵
$\large {S_w}={S_1} + {S_2}$
类间离散度矩阵：
$\large {S_b}=（\mu_1-\mu_2）（\mu_1-\mu_2）^T$

2.一维y空间
各类样本的均值:
$\large \tilde{\mu_i} = \frac{ 1}{N_i } \sum_{y\in \omega_i} y =\frac{ 1}{N_i } \sum_{x\in \omega_i} w^Tx$
$\large \tilde{\mu_i} = w^T \mu_i$
类内离散度
$\large \tilde{S_i}^2=\sum_{y\in \omega_i} (y-\tilde\mu_i)^2$
$\large \tilde{S_i}^2=\sum_{x\in \omega_i} (w^T x-w^T \mu_i)^2$
$\large \tilde{S_i}^2=\sum_{x\in \omega_i} w^T( x-\mu_i)( x-\mu_i)^Tw$
$\large \tilde{S_i}^2=w^TS_iw$
总类内离散度
$\large \tilde{S_w}=\tilde{S_1}^2 + \tilde{S_2}^2$
$\large \tilde{S_w}=w^T(S_1 + S_2)w$
$\large \tilde{S_w}=w^T(S_w)w$

最终：
$\Large J(w) = \frac { w^T S_b w }{w^T S_w w }$

2. 线性判别分析（多类情况）

有C个类别，需要K维向量（或者叫做基向量）来做投影。
$\large W = [w_1 w_2 ...w_k]^T$
$\large y = W^T x$

1.d维x空间
类内离散度矩阵：
$\large {S_i}=\sum_{x\in \omega_i} ( x-\mu_i)( x-\mu_i)^T$
总类内离散度矩阵
$\large {S_w}=\sum_{i=1}^{c} S_i$
类间离散度矩阵(每类均值点相对于样本中心的散列情况)：
$\large {S_b}= \sum_{i=1}^{c}N_i（\mu_i-\mu）（\mu_i-\mu）^T$
其中$\mu$为所有样本的均值：
$\large \mu= \frac{1}{N}\sum_{\forall x}x = \frac{1}{N} \sum_{x\in w_i}N_i \mu_i$

2.k维y空间
$\large \tilde{\mu_i} = \frac{ 1}{N_i } \sum_{y\in \omega_i} y$
$\large \tilde{\mu} = \frac{ 1}{N } \sum_{\forall y} y$
类内离散度矩阵：
$\large \tilde{S_i}=\sum_{y\in \omega_i} ( y-\tilde{\mu_i})( y-\tilde{\mu_i})^T$
总类内离散度
$\large \tilde{S_w}=\sum_{i=1}^{c} \tilde{S_i} = \sum_{i=1}^{c} \sum_{y\in \omega_i} ( y-\tilde{\mu_i})( y-\tilde{\mu_i})^T$
类间离散度矩阵(每类均值点相对于样本中心的散列情况)：
$\large \tilde{S_b}= \sum_{i=1}^{c}N_i（\tilde\mu_i-\tilde\mu）（\tilde\mu_i-\tilde\mu）^T$

最终：
$\Large J(w) = \frac { |\tilde S_b|}{|\tilde S_w |} = \frac { w^T S_b w }{w^T S_w w }$

3. 求解（Fisher linear discrimination）

利用Lagrange乘法算子求解：
$L(w, \lambda) = w^TS_bw - \lambda( w^TS_ww -c)$
令 $\frac{\partial L}{\partial w} = S_bw - \lambda S_w w = 0$
$S_b w^* - \lambda S_w w^* = 0$
即：
$S_b w^* = \lambda S_w w^*$
$S_w^{-1} S_b w^* = \lambda w^*$
$\large \because {S_b}=（\mu_1-\mu_2）（\mu_1-\mu_2）^T$
$\large \therefore {S_b}w^*=（\mu_1-\mu_2）（\mu_1-\mu_2）^T w^* = （\mu_1-\mu_2）\lambda_w$ 其中$\lambda_w$为常数
$\large \therefore S_w^{-1} （\mu_1-\mu_2）\lambda_w = \lambda w^*$
由于对w扩大缩小任何倍不影响结果，因此可以约去两边的未知常数$\lambda_w 和 \lambda$，得到
$\large S_w^{-1} （\mu_1-\mu_2）= w^*$

代码实现：

int cv_lda(void)  
{  
    //sampledata  
    double sampledata[6][2]={{0,1},{0,2},{2,4},{8,0},{8,2},{9,4}};  
    Mat mat=Mat(6,2,CV_64FC1,sampledata);  
    //labels  
    vector<int>labels;  

    for(int i=0;i<mat.rows;i++)  
    {  
        if(i<mat.rows/2)  
        {  
            labels.push_back(0);  
        }  
        else  
        {  
            labels.push_back(1);  
        }  
    }  

    //do LDA  
    //初始化并计算，构造函数中带有计算  
    LDA lda=LDA(mat,labels,1);  
    //get the eigenvector  
    //获得特征向量  
    Mat eivector=lda.eigenvectors().clone();  

    cout<<"特征向量（double）类型:"<<endl;  
    for(int i=0;i<eivector.rows;i++)  
    {  
        for(int j=0;j<eivector.cols;j++)  
        {  
            cout<<eivector.ptr<double>(i)[j]<<" ";  
        }  
        cout<<endl;  
    }  


    //------------------------------计算两个类心------------  
    //针对两类分类问题，计算两个数据集的中心  
    int classNum=2;  
    vector<Mat> classmean(classNum);  
    vector<int> setNum(classNum);  

    for(int i=0;i<classNum;i++)  
    {  
        classmean[i]=Mat::zeros(1,mat.cols,mat.type());  //初始化类中均值为0  
        setNum[i]=0;  //每一类中的条目数  
    }  

    Mat instance;  
    for(int i=0;i<mat.rows;i++)  
    {  
        instance=mat.row(i);//获取第i行  
        if(labels[i]==0)  //如果标签为0  
        {     
            add(classmean[0], instance, classmean[0]);  //矩阵相加  
            setNum[0]++;  //数量相加  
        }  
        else if(labels[i]==1)  //对于第1类的处理  
        {  
            add(classmean[1], instance, classmean[1]);  
            setNum[1]++;  
        }  
        else  
        {}  
    }  
    for(int i=0;i<classNum;i++)   //计算每一类的均值  
    {  
        classmean[i].convertTo(classmean[i],CV_64FC1,1.0/static_cast<double>(setNum[i]));  
    }  
    //----------------------------------END计算类心-------------------------  


    vector<Mat> cluster(classNum);  //一共2类  


    cout<<"特征向量："<<endl;  
    cout<<eivector<<endl;   //此时的特征向量是一个列向量  


    cout<<endl<<endl;  
    cout<<"第一种方式(手动计算)："<<endl;  
    //1.投影的第一种方式：Y=X*W  
    //有的教程写成Y=W^T*X,（此时的X是列向量看待的所以需要将w转置）  
    for(int i=0;i<classNum;i++)  
    {  
        cluster[i]=Mat::zeros(1,1,mat.type()); //初始化0  
        //特征向量的转置同类均值相乘)  
        cluster[i]=classmean[i]*eivector;  
    }  

    cout<<"The project cluster center is:"<<endl;  //计算均值的投影  
    for(int i=0;i<classNum;i++) //输出两类中心的投影值  
    {  
        cout<<cluster[i].at<double>(0)<<endl;  
    }  

    //2.第二种方式使用内置函数计算  
    //第一个中心  
    cout<<endl<<"第二种方式:";  
    cout<<endl<<"第一个类均值的投影:"<<endl;  
    cout<<lda.project(classmean[0]).at<double>(0);  
    cout<<endl<<"第二个类均值的投影"<<endl;  
    cout<<lda.project(classmean[1]).at<double>(0);  


    system("pause");  
    return 0;  
}  

结果：

reference:
[1]. http://www.cnblogs.com/jerrylead/archive/2011/04/21/2024384.html
[2]. http://docs.opencv.org/modules/contrib/doc/facerec/facerec_tutorial.html