2019-04-01

数字媒体技术--求矩阵的对角阵

发现ppt上面的思路有错，结果照着它的方式实现就错了。

不过也有意外的收获，相当于是自己手写了一遍矩阵的对角化的变换。

要变换的矩阵保存在A矩阵中。

然后坐标变换的矩阵保存在C矩阵中。其中变换的原理可以看西电版的线性代数里面的计算。

#include <bits/stdc++.h>
using namespace std;
const int maxn = 100+10;
double eps = 1e-7;
struct Mat{
    double mat[maxn][maxn];
    int row, col;
    Mat(int _row, int _col){
        memset(mat, 0, sizeof(mat));
        for(int i=0; i<maxn; i++) mat[i][i] = 1;
        row = _row, col = _col;
    }
    Mat(){
        memset(mat, 0, sizeof(mat));
        for(int i=0; i<maxn; i++){mat[i][i] = 1;}
    }
};
struct Matrix{
    Mat a;
    Matrix mul(Mat b) {
        Matrix temp;
        memset(temp.a.mat, 0, sizeof(temp.a.mat));
        for(int i=0; i<a.row; i++){
            for(int j=0; j<b.col; j++){
                for(int k=0; k<a.col; k++){
                    temp.a.mat[i][j] += a.mat[i][k]*b.mat[k][j];
                }
            }
        }
        temp.a.row = a.row, temp.a.col = b.col;
        return temp;
    }
    Mat trans(){
        Mat temp = a;
        for(int i=0; i<a.row; i++){
            for(int j=0; j<a.col; j++){
                temp.mat[j][i] = a.mat[i][j];
            }
        }
        swap(a.row, a.col);
        memset(a.mat, 0, sizeof(a.mat));
        for(int i=0; i<a.row; i++){
            for(int j=0; j<a.col; j++){
                a.mat[i][j] = temp.mat[i][j];
            }
        }
    }
    Matrix operator / (double di){
        Matrix temp;
        for(int i=0; i<a.row; i++){
            for(int j=0; j<a.col; j++){
                temp.a.mat[i][j] = a.mat[i][j]/di;
            }
        }
        temp.a.row = a.row, temp.a.col = a.col;
        return temp;
    }
};
void prin(Matrix a){
    for(int i=0; i<a.a.row; i++){
        for(int j=0; j<a.a.col; j++){
            printf("%lf ", a.a.mat[i][j]);
        }
        printf("\n");
    }
}
//double table[maxn][maxn] =
//            {{10, 15, 23, 11, 42, 9, 11, 8, 11, 21},
//            {15, 46, 21, 9, 45, 48, 21, 5, 12, 20},
//            {29, 13, 30, 35, 11, 5, 14, 15, 21, 25}};
double table[maxn][maxn] = {{35, 1, 13, 2, 4, 41, 34, 15, 47, 1},
                            {21, 19, 38, 39, 9, 24, 22, 32, 35, 37},
                            {13, 33, 32, 8, 5, 24, 47, 17, 29, 11}};
//double table[maxn][maxn] = {{31.4362,-8.8665, 3.9004, 0.4302, 2.6514,34.0677,19.8268,10.4482,38.7669,-1.5558},
//{20.6173,18.0285,37.0580,38.7645, 8.8528,23.2935,20.6164,31.4995,34.1463,36.6762},
//{13.0000,33.0000,32.0000, 8.0000, 5.0000,24.0000,47.0000,17.0000,29.0000,11.0000}};
int row, col;
double A[maxn][maxn];
double C[maxn][maxn];
void up_train(int row, int col){
    memset(C, 0, sizeof(C));
    cout<<endl<<endl<<endl<<endl;
    cout<<row<<" "<<col<<endl;
    for(int i=0; i<max(row, col); i++) C[i][i] = 1;
    for(int i=0; i<row; i++){
        if(fabs(A[i][i])<eps&&i==0){
            int k=i+1;
            while(fabs(A[k][i])<eps&&k<row){
                k++;
            }
            bool has = false;
            if(fabs(A[k][i])>eps&&k<row) has = true;
            //cout<<"in  "<<k<<" "<<A[k][i]<<endl;
            if(has){
                for(int t=0; t<col; t++) {
                    A[i][t] += A[k][t];
                    //C[i][t] += C[k][t];
                }
                for(int t=0; t<row; t++){
                    A[t][i] += A[t][k];
                    //C[t][i] += C[t][k];
                }
                for(int t=0; t<row; t++){
                    //A[t][i] += A[t][k];
                    C[t][i] += C[t][k];
                }
            }
            else{
                for(int t=0; t<col; t++) {
                    A[i][t] += A[i-1][t];
                    //C[i][t] += C[i-1][t];
                }
                for(int t=0; t<col; t++) {
                    A[t][i] += A[t][i-1];
                    //C[t][i] += C[t][i-1];
                }
                for(int t=0; t<col; t++) {
                    //A[t][i] += A[t][i-1];
                    C[t][i] += C[t][i-1];
                }
            }
        }
        for(int j=0; j<i; j++){
            double tim = A[i][j]/A[j][j];
            A[i][j] = 0.0;
            //C[i][j] -= tim*C[j][j];
            for(int k=j+1; k<col; k++){
                A[i][k] -= tim*A[j][k];
                //C[i][k] -= tim*C[j][k];
            }
            for(int k=0; k<row; k++){
                A[k][i] -= tim*A[k][j];
                //C[k][i] -= tim*C[k][j];
            }
            for(int k=0; k<row; k++){
                //A[k][i] -= tim*A[k][j];
                C[k][i] -= tim*C[k][j];
            }
            //cout<<i<<",--"<<j<<" "<<A[i][j]<<endl;
        }
        for(int i=0; i<row; i++){
            for(int j=0; j<col; j++){
                cout<<A[i][j]<<" ";
            }
            cout<<endl;
        }
    }
    printf("-----------\n");
    for(int i=0; i<row; i++){
        for(int j=0; j<col; j++){
            printf("%lf ", C[i][j]);
        }
        printf("\n");
    }
    for(int i=row-1; i>=0; i--){
        for(int j=col-1; j>i; j--){
            if(A[j][j] == 0) continue;
            //cout<<"kuku "<<i<<" "<<j<<endl;
            double tim = A[i][j]/A[j][j];
            A[i][j] = 0.0;
            for(int k=j-1; k>i; k--){
                A[i][k] -= tim*A[j][k];
                //C[i][k] -= tim*C[j][k];
            }
            for(int k=0; k<col; k++){
                A[k][i] -= tim*A[k][j];
                //C[k][i] -= tim*C[k][j];
            }
            for(int k=0; k<col; k++){
                //A[k][i] -= tim*A[k][j];
                C[k][i] -= tim*C[k][j];
            }
            //cout<<i<<",--"<<j<<" "<<A[i][j]<<endl;
        }
    }
    printf("second time:---------\n");
    for(int i=0; i<row; i++){
        for(int j=0; j<col; j++){
            printf("%lf ", C[i][j]);
        }
        printf("\n");
    }
    printf("seperate line\n");
    for(int i=0; i<row; i++){
        for(int j=0; j<col; j++){
            printf("%lf ", A[i][j]);
        }
        printf("\n");
    }
}
double te[maxn][maxn] = {{0, 1.0/2, -5.0/2}, {1.0/2, 0, 1.0/2}, {-5.0/2, 1.0/2, 0}};
void test(){
    row = col = 3;
    for(int i=0; i<row; i++){
        for(int j=0; j<col; j++){
            A[i][j] = te[i][j];
        }
    }
    for(int i=0; i<row; i++){
        for(int j=0; j<col; j++){
            cout<<A[i][j]<<" ";
        }
        cout<<endl;
    }
}
int main(){
    ios::sync_with_stdio(false);
    Matrix a1;
    a1.a.row = 3, a1.a.col = 10;
    for(int i=0; i<a1.a.row; i++){
        for(int j=0; j<a1.a.col; j++){
            a1.a.mat[i][j] = table[i][j];
        }
    }
    for(int i=0; i<a1.a.row; i++){
        double ave = 0;
        for(int j=0; j<a1.a.col; j++){
            ave += a1.a.mat[i][j];
        }
        cout<<ave<<" ";
        ave /= a1.a.col;
        cout<<ave<<endl;
        for(int j=0; j<a1.a.col; j++){
            a1.a.mat[i][j] -= ave;
        }
    }
    Matrix a_t = a1;
    a_t.trans();
    prin(a_t);
    Matrix ans = (a1.mul(a_t.a));
    //求出协方差矩阵，注意除的是n-1;
    ans = ans/(a1.a.col-1);
    printf("\n\n\n");
    prin(ans);
    row = ans.a.row;
    col = ans.a.col;
    for(int i=0; i<row; i++){
        for(int j=0; j<col; j++){
            A[i][j] = ans.a.mat[i][j];
        }
    }
    //test();
    up_train(row, col);
    return 0;
}

pca为什么取特征值最大的列证明

link

本帖最后由 miniwhale 于 2011-1-14 20:39 编辑
本科数学的问题之一就是没有把线性代数放在应有的地位上，没有讲透。而大家都知道后续科目中的公式基本上都要靠矩阵表达，这给后续学习带来不少困难，当初我也是重新自学的。闲话少叙，进入正题
绝大多数“多元统计”教材讲到多元正态分布时，都有这么一道习题，包括相关的2个证明：
1)、Y是随机变量，X是随机变量，c是常量，设Y=c.X，证明D(Y)=c^2*D(X)
2)、Y是随机变量，x是随机向量，c是常向量，设Y=c’x，证明D(Y)=c’.D(x).c
利用此结果：
∵第i个主成分Yi=Vi’.X,这里Vi是D(x)的第i个特征向量
∴D(Yi)=Vi’.D(X).Vi
∴D(Yi)=Vi’.[D(X).Vi]=Vi’[λi.Vi]，第2个等号的依据是特征向量的定义
∴D(Yi)=λi[Vi’.Vi]=λi，第1个等号的依据是提取常系数，第2个等号的依据是Vi是单位向量
All done。
ps：
其实，如果完全用线性代数的语言，完全可以把主成分分析表达的更漂亮
上面的习题有时还会有第3个证明：3)、y是随机向量，x是随机向量，C是常矩阵，设y=C’x，证明D(y)=C’.D(x).C
利用此结果：
主成分y=V’.x，这里V是特征向量排列成的矩阵
∴D(y)=V’.D(x).V=Λ，这里Λ是特征值沿着对角线排列成的矩阵

未解决的问题

本文标题:数字媒体技术--求矩阵的对角阵

文章作者:Babydragon

发布时间:2019-04-01, 19:32:20

最后更新:2019-04-01, 23:25:01

原始链接:http://baolintian.github.io/2019/04/01/数字媒体技术-求矩阵的对角阵/

许可协议: "署名-非商用-相同方式共享 4.0" 转载请保留原文链接及作者。