精华内容
下载资源
问答
  • 2022-03-13 22:12:28

    C语言实现(复数)矩阵运算

    代码并不是很严谨,没有考虑内存释放之类的问题,初次尝试还有待完善,提供矩阵运算思路、希望能够帮助到有需要的同学。

    1.首先创建复数类型结构体

    typedef struct
    {
        float real;
        float imag;
    }Complex_F;
    

    结构体Complex_F中定义了复数的实部与虚部,这里我使用Complex_F *A或者Complex_F A[row*col]来表示一个矩阵,由于是行排列,所以没有使用二维数组,只需知道A[i][j]=A[i*col+j]即可。
    下面是所使用的代码,重点是LU分解求逆矩阵,
    参考了https://blog.csdn.net/weixin_41645983/

    代码实现功能

    代码实现了复数的加减乘除,复数矩阵的初始化、拷贝、打印,复数矩阵的乘法、求逆以及转置(复数矩阵一般用共轭转置情形较多,比如要求H矩阵每列范数,可用H共轭转置乘以H,取对角线元素即可)。

    #include<stdio.h>
    #include<math.h>
    #include<stdlib.h>
    typedef struct
    {
        float real;
        float imag;
    }Complex_F;//我们定义一个复数类型,构建Complex_F类型数组
    
    void initMatrix(Complex_F* A, int row, int col);
    Complex_F MulComplex(const Complex_F a, const Complex_F b);
    Complex_F addComplex(const Complex_F a, const Complex_F b);
    void addMatrix(const Complex_F* A, const Complex_F* B, Complex_F* C, int row, int col);
    void PrintMatrix(const Complex_F* A, int row, int col);
    void MulMatrix(const Complex_F* A, const Complex_F* B, Complex_F* C, int row, int AB_rc, int col);
    void TransMatrix(const Complex_F* A, Complex_F* B, int row, int col);
    void MatrixCopy(const Complex_F* A, Complex_F* B, int row, int col);
    void Inverse_LU(Complex_F* A, int row, int col);
    
    
    void initMatrix(Complex_F* A, int row, int col)
    {
        int i, j;
        for (i = 0; i < row; i++)
        {
            for (j = 0; j < col; j++)
            {
                A[i * col + j].real = 0.0;
                A[i * col + j].imag = 0.0;
            }
        }
        //以上写法是将数组看成矩阵,行排列读取第i行第j列用i*col+j;
    }
    //当传进去的矩阵元素在函数中不发生改变时,用const Complex_F* A
    //只要不想改变实参时,都可以加上const避免出错
    Complex_F MulComplex(const Complex_F a, const Complex_F b)
    {
        Complex_F c;
        c.real = a.real * b.real - a.imag * b.imag;
        c.imag = a.real * b.imag + a.imag * b.real;
        return c;
    }
    Complex_F addComplex(const Complex_F a, const Complex_F b)
    {
        Complex_F c;
        c.real = a.real + b.real;
        c.imag = a.imag + b.imag;
        return c;
    }
    Complex_F subComplex(const Complex_F a, const Complex_F b)
    {
        Complex_F c;
        c.real = a.real - b.real;
        c.imag = a.imag - b.imag;
        return c;
    }
    Complex_F divComplex(const Complex_F a, const Complex_F b)
    {
        Complex_F c;
        float norm = b.real * b.real + b.imag * b.imag;
        c.real = (a.real * b.real + a.imag * b.imag) / norm;
        c.imag = (a.imag * b.real - a.real * b.imag) / norm;
        return c;
    }
    void addMatrix(const Complex_F* A, const Complex_F* B, Complex_F* C, int row, int col)
    //传递指针可以不用带返回值,我们想得到C的结果,在传入进去前给C分配空间即可,A和B加上const
    {
        int i, j;
        for (i = 0; i < row; i++) {
            for (j = 0; j < col; j++)
            {
                C[i * col + j].real = A[i * col + j].real + B[i * col + j].real;
                C[i * col + j].imag = A[i * col + j].imag + B[i * col + j].imag;
            }
        }
    }
    void PrintMatrix(const Complex_F* A, int row, int col)
    {
        int i, j;
        for (i = 0; i < row; i++)
        {
            for (j = 0; j < col; j++)
            {
                printf("%.3f+(%.3f)i\t", A[i * col + j].real, A[i * col + j].imag);
            }
            printf("\n");
        }
    }
    
    void MulMatrix(const Complex_F* A, const Complex_F* B, Complex_F* C, int row, int AB_rc, int col)
    {
        int i, j, k;
        Complex_F sum, tmp;
        for (i = 0; i < row; i++)
        {
            for (j = 0; j < col; j++)
            {
    
                sum.real = 0.0;
                sum.imag = 0.0;
                for (k = 0; k < AB_rc; k++)
                {
                    sum = addComplex(sum, MulComplex(A[i * AB_rc + k], B[k * col + j]));
                }
                C[i * col + j].real = sum.real;
                C[i * col + j].imag = sum.imag;
            }
        }
    }
    
    void MatrixCopy(const Complex_F* A, Complex_F* B, int row, int col)
    {
        int i, j;
        for (i = 0; i < row; i++)
        {
            for (j = 0; j < col; j++)
            {
                B[i * col + j] = A[i * col + j];
            }
        }
    }
    void TransMatrix(const Complex_F* A, Complex_F* B, int row, int col)
    {
        int i, j;
        for (i = 0; i < row; i++)
        {
            for (j = 0; j < col; j++)
            {
                B[j * row + i] = A[i * col + j];
            }
        }
    }
    void TransConjMatrix(const Complex_F* A, Complex_F* B, int row, int col)
    {
        int i, j;
        for (i = 0; i < row; i++)
        {
            for (j = 0; j < col; j++)
            {
                B[j * row + i].real = A[i * col + j].real;
                B[j * row + i].imag = -A[i * col + j].imag;
            }
        }
    }
    //LU分解求逆矩阵
    //LU分解将一个矩阵分解成一个单位下三角矩阵和一个上三角矩阵的乘积
    //我们希望L的对角线上的元素为1,使用Doolittle分解得到
    //对于满秩矩阵A可以通过相乘一个消元矩阵得到一个上三角矩阵U
    //A=LU div(A)=div(U)div(L)
    
    void Inverse_LU(Complex_F* A, int row, int col)
    {
        Complex_F* L;
        Complex_F* U;
        Complex_F* L_Inverse;
        Complex_F* U_Inverse;
        L = (Complex_F*)malloc(row * col * sizeof(Complex_F));
        U = (Complex_F*)malloc(row * col * sizeof(Complex_F));
        L_Inverse = (Complex_F*)malloc(row * col * sizeof(Complex_F));
        U_Inverse = (Complex_F*)malloc(row * col * sizeof(Complex_F));
        initMatrix(L, row, col);
        initMatrix(U, row, col);
        initMatrix(L_Inverse, row, col);
        initMatrix(U_Inverse, row, col);
        int i, j, k, t;
        Complex_F s;
        for (i = 0; i < row; i++)//对角元素
        {
            L[i * col + i].real = 1.0;
            L[i * col + i].imag = 0.0;
        }
        for (j = 0; j < col; j++)
        {
            U[j] = A[j];
        }
        for (i = 1; i < col; i++)
        {
            L[i * col] = divComplex(A[i * col], U[0]);
    
        }
        for (k = 1; k < row; k++)
        {
            for (j = k; j < col; j++)
            {
                s.imag = 0.0;
                s.real = 0.0;
                for (t = 0; t < k; t++) {
                    s = addComplex(s, MulComplex(L[k * col + t], U[t * col + j]));
                }
                U[k * col + j] = subComplex(A[k * col + j], s);
            }
            for (i = k; i < col; i++)
            {
                s.imag = 0.0;
                s.real = 0.0;
                for (t = 0; t < k; t++)
                {
                    s = addComplex(s, MulComplex(L[i * col + t], U[t * col + k]));
                }
                L[i * col + k] = divComplex(subComplex(A[i * col + k], s), U[k * col + k]);
            }
        }
    
        for (i = 0; i < row; i++)
        {
            L_Inverse[i * col + i].imag = 0.0;
            L_Inverse[i * col + i].real = 1.0;
        }
        for (j = 0; j < col; j++)
        {
            for (i = j + 1; i < row; i++)
            {
                s.imag = 0.0;
                s.real = 0.0;
                for (k = j; k < i; k++) {
                    s = addComplex(s, MulComplex(L[i * col + k], L_Inverse[k * col + j]));
                }
                L_Inverse[i * col + j].real = (-1) * MulComplex(L_Inverse[j * col + j], s).real;
                L_Inverse[i * col + j].imag = (-1) * MulComplex(L_Inverse[j * col + j], s).imag;
            }
        }
        Complex_F di;
        di.real = 1.0;
        di.imag = 0.0;
        for (i = 0; i < col; i++)                    //按列序,列内按照从下到上,计算u的逆矩阵
        {
            U_Inverse[i * col + i] = divComplex(di, U[i * col + i]);
    
        }
        for (j = 0; j < col; j++)
        {
            for (i = j - 1; i >= 0; i--)
            {
                s.imag = 0.0;
                s.real = 0.0;
                for (k = i + 1; k <= j; k++)
                {
                    s = addComplex(s, MulComplex(U[i * col + k], U_Inverse[k * col + j]));
    
                }
                s.imag = -s.imag;
                s.real = -s.real;
                U_Inverse[i * col + j] = divComplex(s, U[i * col + i]);
            }
        }
        MulMatrix(U_Inverse, L_Inverse, A, row, row, col);
    }
    /*Complex_F Det(const Complex_F* A,int row,int col)
    {
        if(row!=col)
        {
            printf("error:please input a square matrix!\n");
            return;
        }
        else{
    
    
        }
    }*/
    int main()
    {
        Complex_F A[4] = { {1,0},{10,0},{5,0},{2,0} };
        Complex_F B[4] = { {1.0,2.0},{1.0,2.0},{1.1,1.1},{1.0,2.0} };
        Complex_F* D;
        Complex_F* C;
        Complex_F* A_;
        C = (Complex_F*)malloc(4 * sizeof(Complex_F));//空指针一定要记得分配空间!!
        D = (Complex_F*)malloc(4 * sizeof(Complex_F));
        A_ = (Complex_F*)malloc(4 * sizeof(Complex_F));
        //MulMatrix(A,B,C,2,2,2);
        //PrintMatrix(C,2,2);
        //TransConjMatrix(A,D,2,2);//计算A矩阵每列的范数平方可以用A的转置乘以A,取对角元素得到,复数矩阵则是共轭转置
        //MulMatrix(D,A,C,2,2,2);
        //PrintMatrix(C,2,2);
        Inverse_LU(A, 2, 2);
        PrintMatrix(A, 2, 2);
    
        return 0;
    }
    
    更多相关内容
  • 本章开始,笔者将详细讲解矩阵的QR分解、特征值、特征向量等运算,并给出函数和测试demo的C代码,以及与matlab计算的结果; 并且,本章相当于前面几章的大杂烩,前面所有的结构体、宏定义、函数本章基本全部都有用到...

    一、前言

    1. 本章开始,笔者将详细讲解矩阵的QR分解特征值特征向量等运算,并给出函数和测试demo的C代码,以及与matlab计算的结果;
    2. 并且,本章相当于前面几章的大杂烩,前面所有的结构体、宏定义、函数本章基本全部都有用到。
    3. 此外,由于本章求解QR矩阵/特征值/特征向量,实数和复数的过程原理完全一样,读者可以自行修改参数,简化复数计算的部分即可,或者完全不用修改,直接调用我的函数,虚部设置为0退化为实数运算即可;
    4. 另,这将也是我此连载文章的系列精品之一,矩阵运算的功能越来越复杂,但是笔者全部分模块进行,化繁为简,化整为零,所以需要调用的小函数越来越多,读者遇到不熟悉的函数可以随时Jump到(一)(二)(三)(四)的内容。
      (PS:另外这两个编辑器为啥功能不能优势互补,而且都没有贴matlab code的编辑器?)
      在这里插入图片描述

    二、QR分解

    QR分解的数学原理读者可以直接点击链接复习一下,QR分解C语言实现我参考的是此大佬的。笔者这里整理出来了,C code直接给出:

    /* QR Decompose */
    void QR(const Matrix* A, Matrix* Q, Matrix* R)
    {
    	if (IsNullComplexMatrix(A))
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return;
    	}
    	// Not A Square Matrix
    	if (A->row != A->column)
    	{
    		printf("ERROE: Not a square matrix!\n");
    		return;
    	}
    	
    	int i, j, k, m;
    	int size;
    	const int N = MatrixRow(A);
    	ComplexType temp;
    
    	// Column Vector Saving Column Vectors of A
    	Matrix a, b;
    	InitComplexMatrix(&a, N, 1);
    	InitComplexMatrix(&b, N, 1);
    	size = MatrixSize(A);
    	if (MatrixSize(Q) != size)
    	{
    		// free(Q->arrayComplex);
    		DestroyComplexMatrix(Q);                   // Free The Initial Matrix
    		InitComplexMatrix(Q, A->row, A->column);   // Reset Size and Initialize
    	}
    
    	if (MatrixSize(R) != size)
    	{
    		// free(R->arrayComplex);
    		DestroyComplexMatrix(R);
    		InitComplexMatrix(R, A->row, A->column);
    	}
    
    	for (j = 0; j < N; ++j)
    	{
    		for (i = 0; i < N; ++i)
    		{
    			a.arrayComplex[i] = b.arrayComplex[i] = A->arrayComplex[i * A->column + j];   // Cols Vector of A
    		}
    
    		for (k = 0; k < j; ++k)
    		{
    			R->arrayComplex[k * R->column + j]._Val[0] = 0;
    			R->arrayComplex[k * R->column + j]._Val[1] = 0;
    //			InitComplex(R->arrayComplex + k * R->column + j);  // reset
    			for (m = 0; m < N; ++m)
    			{
    				R->arrayComplex[k * R->column + j] = AddComplex(R->arrayComplex[k * R->column + j], \
    					_Cmulcc(a.arrayComplex[m], Q->arrayComplex[m * Q->column + k]));
    			}
    
    			for (m = 0; m < N; ++m)
    			{
    				b.arrayComplex[m] = SubComplex(b.arrayComplex[m], _Cmulcc(R->arrayComplex[k * R->column + j], \
    					Q->arrayComplex[m * Q->column + k]));
    			}
    		}
    
    		temp = MatrixNorm2(&b);
    		R->arrayComplex[j * R->column + j] = temp;
    
    		for (i = 0; i < N; ++i)
    		{
    			Q->arrayComplex[i * Q->column + j] = DivComplex(b.arrayComplex[i], temp);
    		}
    	}
    	// Free Local Memory
    	Matrix ComplexMatrixArray[] = { a, b };
    	int numDoubleMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
    	DestroyComplexMatrixArray(ComplexMatrixArray, numDoubleMatrixArray);      // Complpex Matrix
        // OR
    //	DestroyComplexMatrix(&a);
    //	DestroyComplexMatrix(&b);
    }
    

    说明
    (1) 因为在后面的特征值的求解函数中需要用到,单单QR分解本身数学意义不大,因此单独的测试就免了;
    (2) 后面释放结构体指针(也就是矩阵内存)的函数,调用的是DestroyComplexMatrixArray()而非DestroyComplexMatrix(),前面第一章给了二者实现的代码,读者留意一下,我个人建议读者定义的矩阵不多,最好不要调用这个destroy函数,老老实实一个一个调用DestroyComplexMatrix()/DestroyDoubleMatrix()函数释放内存;
    (3) 另外,注意矩阵的个数numDoubleMatrixArray必须在调用DestroyComplexMatrixArray()函数之前求出来赋值给变量然后传给形参,而不能在DestroyComplexMatrixArray()函数里面求解,且ComplexMatrixArray必须定义为数组而不能是指针,因为无论是指针还是数组,实参传给形参时候,系统都会在函数内部将其当作指针处理,sizeof(ptr) ≡ 4Byte(64bit系统下),而不是指向的那段内存的实际大小,但是sizeof(数组名)却是实际的那段内存的大小,这也是数组和指针的区别,C语言里面这一块也是一个比较容易出错的地方(笔者此处也差点搞忘了!!!)。

    三、矩阵特征值

    特征值计算比较复杂,关于特征值和特征向量的求解,线性代数理论知识大家可以参考此篇文章,关于特征值和特征向量C语言的实现,笔者是根据此位大佬修改的,这里我直接贴出笔者魔改后的代码:

    /* eigen values */
    void EigenValue(const Matrix* matrix, Matrix* eigenvalue)
    {
    	const int NUM = 100;   // Iteration Times
    	
    	// Local Matrice
    	Matrix Q, R, temp;
    	// Initiate
    	InitComplexMatrix(&Q, matrix->row, matrix->column);
    	InitComplexMatrix(&R, matrix->row, matrix->column);
    	InitComplexMatrix(&temp, matrix->row, matrix->column);
    
    	if (IsNullComplexMatrix(matrix))
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return;
    	}
    	// Copy matrix to temp in order not to change matrix
    	CopyMatrix(matrix, &temp);
    	// QR Decompose and 
    	for (int k = 0; k < NUM; ++k)
    	{
    		QR(&temp, &Q, &R);
    		MatrixMulMatrix(&R, &Q, &temp);
    	}
    	// Abstract Eigen Values from the Diagonal Elements of temp =  Q * R
    	for (int k = 0; k < temp.row; ++k)
    	{
    		eigenvalue->arrayComplex[k] = temp.arrayComplex[k * temp.column + k];
    	}
    	// Free Local Memory
    	Matrix ComplexMatrixArray[] = { R, Q, temp };
    	int numComplexMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
    	DestroyComplexMatrixArray(ComplexMatrixArray, numComplexMatrixArray);      // Complpex Matrix
    //	return eigenvalue;
    }
    

    说明
    (1) 因为在后面的特征向量的求解中需要用到,因此单独的测试就免了;
    (2) 同样这里销毁结构体内部指针(即矩阵)内存调用的也是一次销毁函数DestroyComplexMatrixArray();
    (3) 读者注意一下笔者在此特征值求解函数中调用的一些函数,如果忘记了可以跳到前面几章节review。

    三、矩阵特征向量

    我是一个没有感情的复读机:特征向量的计算更为复杂,关于特征向量的求解,线性代数理论知识大家可以参考此篇文章,关于特征向量C语言的实现,笔者是根据此位大佬修改的,这里我直接贴出笔者修魔改后的代码:

    /* Negative Complex : Complex_B = -Complex_A = -creal(Complex_A) - cimag(Complpex_A) */
    ComplexType NegativeComplex(const ComplexType Complex_A)
    {
    	ComplexType Complex_B;
    	Complex_B = _Cmulcr(Complex_A, -1.0);
    	// OR
    /*
    	Complex_B._Val[0] = -creal(Complex_A);
    	Complex_B._Val[1] = cimag(Complex_A) * (-1.0);
    */
    	return Complex_B;
    }
    
    /* eigen vectors */
    void EigenVector(const Matrix* matrix, const Matrix* eigenvalue, Matrix* eigenvector)
    {
    	if (IsNullComplexMatrix(matrix) || IsNullComplexMatrix(eigenvalue))
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return;
    	}
    	int i, j, q;
    	int m;
    
    	// Access to Eigen Values
    	int count;
    	int num = MatrixRow(matrix);   // = matrix->row: Numbers of Eigen Values or Cols
    	ComplexType evalue;
    
    	// Access to temp
    	ComplexType sum, midsum, mid;
    	Matrix temp;   // temp = A - λI: (A - λI) * x = 0
    	InitComplexMatrix(&temp, matrix->row, matrix->column);
    
    	for (count = 0; count < num; ++count)
    	{
    		// Calculate x: Ax = λ * x
    		evalue = eigenvalue->arrayComplex[count];
    		CopyMatrix(matrix, &temp);
    		for (i = 0; i < temp.column; ++i)
    		{
    			temp.arrayComplex[i * temp.column + i] = SubComplex(temp.arrayComplex[i * temp.column + i], evalue);
    			//			temp->arrayComplex[i * temp->column + i]._Val[0] -= creal(evalue);
    			//			temp->arrayComplex[i * temp->column + i]._Val[1] -= cimag(evalue);
    		}
    
    		// Transform temp to Ladder Matrix
    		for (i = 0; i < temp.row - 1; ++i)
    		{
    			mid._Val[0] = creal(temp.arrayComplex[i * temp.column + i]);   // Diagonal Element
    			mid._Val[1] = cimag(temp.arrayComplex[i * temp.column + i]);
    			for (j = i; j < temp.column; ++j)
    			{
    				temp.arrayComplex[i * temp.column + j] = DivComplex(temp.arrayComplex[i * temp.column + j], mid);
    			}
    
    			for (j = i + 1; j < temp.row; ++j)
    			{
    				mid = temp.arrayComplex[j * temp.column + i];
    				for (q = i; q < temp.column; ++q)
    				{
    					temp.arrayComplex[j * temp.column + q] = SubComplex(temp.arrayComplex[j * temp.column + q], \
    						_Cmulcc(temp.arrayComplex[i * temp.column + q], mid));
    
    				}
    			}
    		}
    		midsum._Val[0] = 1; 
    		midsum._Val[1] = 0;
    		eigenvector->arrayComplex[(eigenvector->row - 1) * eigenvector->column + count]._Val[0] = 1;
    		eigenvector->arrayComplex[(eigenvector->row - 1) * eigenvector->column + count]._Val[1] = 0;
    		for (m = temp.row - 2; m >= 0; --m)
    		{
    //			InitComplex(&sum);
    			sum._Val[0] = 0; sum._Val[1] = 0;   // Zero Complex
    			for (j = m + 1; j < temp.column; ++j)
    			{
    				sum = AddComplex(sum, _Cmulcc(temp.arrayComplex[m * temp.column + j],
    					eigenvector->arrayComplex[j * eigenvector->column + count]));
    			}
    			sum = DivComplex(NegativeComplex(sum), temp.arrayComplex[m * temp.column + m]);  // Warning: Parameters' Type
    			//sum = -sum / *(temp.arrayComplex[m * temp.column + m]);
    			midsum = AddComplex(midsum, _Cmulcc(sum, sum));
    			eigenvector->arrayComplex[m * eigenvector->column + count] = sum;
    		}
    
    		midsum = csqrt(midsum);
    		for (i = 0; i < eigenvector->row; ++i)  // One Column Vector--Eigen Vector
    		{
    			eigenvector->arrayComplex[i * eigenvector->column + count] = 
    				DivComplex(eigenvector->arrayComplex[i * eigenvector->column + count], midsum);
    		}
    	}
    	DestroyComplexMatrix(&temp);
    //	return eigenvector;
    }
    

    说明
    (1) 第一个函数用于求解复数的相反数,即a+bi—> -a-bi,当然也可以直接调用complex.h头文件里的_Cmulcr()函数:

    sum = _Cmulcr(sum, -1.0);
    

    (2) 这里只有一个local matrix variable,因此笔者直接调用的DestroyComplexMatrix()函数销毁内存;
    (3) 笔者注释一定要看,包括前面的特征值计算,这样对比数学原理,函数理解起来非常简单。
    (4) 将本章的QR分解、特征值、特征向量函数放在一起测试,测试的demo如下:

    #include<stdio.h>
    #include<stdlib.h>
    #include<complex.h>
    #include <math.h>
    
    typedef _Dcomplex ComplexType;
    typedef double DoubleType;
    
    typedef struct
    {
    	int row, column;
    	ComplexType* arrayComplex;
    }Matrix;
    
    typedef struct
    {
    	int row, column;
    	DoubleType* arrayDouble;
    }Matrix2Double;
    
    typedef enum
    {
    	False = 0, True = 1
    }Bool;
    
    /*
      Initiation of Complex Number
     */
    void InitComplex(ComplexType* Complex)
    {
    	Complex->_Val[0] = 0.0;
    	Complex->_Val[1] = 0.0;
    
    }
    
    /*
      Validity of Complex Matrix
    */
    Bool IsNullComplexMatrix(const Matrix* matrix)
    {
    	int size = matrix->row * matrix->column;
    
    	if (size <= 0 || matrix->arrayComplex == NULL)
    	{
    		return True;
    	}
    	return False;
    }
    
    
    /*
      Initiation of Complex Matrix
     */
    void InitComplexMatrix(Matrix* matrix, int row, int column)
    {
    	int size = row * column * sizeof(ComplexType);
    	if (size <= 0)
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return;
    	}
    	matrix->arrayComplex = (ComplexType*)malloc(size);
    	if (matrix->arrayComplex)
    	{
    		matrix->row = row;
    		matrix->column = column;
    		for (int row_i = 0; row_i < row; ++row_i)
    		{
    			for (int column_j = 0; column_j < column; ++column_j)
    			{
    				InitComplex(matrix->arrayComplex + row_i * matrix->column + column_j);
    
    			}
    		}
    
    	}
    }
    
    /*
      Free Memory of Complex Matrix
    */
    void DestroyComplexMatrix(Matrix* matrix)
    {
    	if (!IsNullComplexMatrix(matrix))
    	{
    		free(matrix->arrayComplex);
    		matrix->arrayComplex = NULL;
    	}
    	matrix->row = matrix->column = 0;
    }
    
    /*
      Free Memory of Complex Matrice Array
    */
    void DestroyComplexMatrixArray(Matrix matrixArray[], int num)  // Array Transfer--->Pointer Transfer
    {
    	if (num)      // if no cell
    	{
    		for (int i = 0; i < num; i++)
    		{
    			DestroyComplexMatrix(&matrixArray[i]);  // Nested Call of DestroyComplexMatrix()
    		}
    	}
    }
    
    /*
      return matrix row size
    */
    int MatrixRow(const Matrix* matrix)
    {
    	return matrix->row;
    }
    
    /*
      return matrix column size
    */
    int MatrixColumn(const Matrix* matrix)
    {
    	return matrix->column;
    }
    
    /*
      return matrix column size
    */
    int MatrixSize(const Matrix* matrix)
    {
    	return (matrix->row * matrix->column);
    }
    
    /*
       Add Complex: Complex_C = Complex_A + Complex_B
    */
    ComplexType AddComplex(const ComplexType Complex_A, const ComplexType Complex_B)
    {
    	ComplexType Complex_C;
    	Complex_C._Val[0] = creal(Complex_A) + creal(Complex_B);
    	Complex_C._Val[1] = cimag(Complex_A) + cimag(Complex_B);
    	return Complex_C;
    }
    
    /*
       Subvision Complex: Complex_C = Complex_A - Complex_B
    */
    ComplexType SubComplex(const ComplexType Complex_A, const ComplexType Complex_B)
    {
    	ComplexType Complex_C;
    	Complex_C._Val[0] = creal(Complex_A) - creal(Complex_B);
    	Complex_C._Val[1] = cimag(Complex_A) - cimag(Complex_B);
    	return Complex_C;
    }
    
    /*
       Division Complex : Complex_C = Complex_A / Complex_B
    */
    ComplexType DivComplex(const ComplexType Complex_A, const ComplexType Complex_B)   // Tip: Return Value Is NOT A Pointer
    {
    	ComplexType Complex_C;
    	Complex_C._Val[0] = (creal(Complex_A) * creal(Complex_B) + cimag(Complex_A)
    		* cimag(Complex_B)) / (pow(creal(Complex_B), 2) + pow(cimag(Complex_B), 2));
    	Complex_C._Val[1] = (cimag(Complex_A) * creal(Complex_B) - creal(Complex_A)
    		* cimag(Complex_B)) / (pow(creal(Complex_B), 2) + pow(cimag(Complex_B), 2));
    	return Complex_C;
    }
    
    /*
       Negative Complex : Complex_B = -Complex_A = -creal(Complex_A) - cimag(Complpex_A)
    */
    ComplexType NegativeComplex(const ComplexType Complex_A)
    {
    	ComplexType Complex_B;
    	Complex_B = _Cmulcr(Complex_A, -1.0);
    	return Complex_B;
    }
    
    /*
       2-norm of a Matrix
    */
    ComplexType MatrixNorm2(const Matrix* matrix)
    {
    	ComplexType norm;
    	norm._Val[0] = 0; norm._Val[1] = 0;
    	if (IsNullComplexMatrix(matrix))
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return (csqrt(norm));
    	}
    	else
    	{
    		for (int i = 0; i < matrix->row; ++i)
    		{
    			for (int j = 0; j < matrix->column; ++j)
    				norm = AddComplex(norm, _Cmulcc(matrix->arrayComplex[i * matrix->column + j], matrix->arrayComplex[i * matrix->column + j]));
    		}
    		return (csqrt(norm));
    	}
    }
    
    /*
      Copy: Complex matrixA = Complex matrixB
    */
    void CopyMatrix(const Matrix* matrix_A, Matrix* matrix_B)
    {
    	if (IsNullComplexMatrix(matrix_A))
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return;
    	}
    	else
    	{
    		int index = 0;
    		for (int row_i = 0; row_i < matrix_A->row; row_i++)
    		{
    			for (int column_j = 0; column_j < matrix_A->column; column_j++)
    			{
    				index = matrix_B->column * row_i + column_j;
    				matrix_B->arrayComplex[index] = matrix_A->arrayComplex[index];
    			}
    		}
    	}
    }
    
    /*
       QR Decompose
    */
    void QR(const Matrix* A, Matrix* Q, Matrix* R)
    {
    	if (IsNullComplexMatrix(A))
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return;
    	}
    	if (A->row != A->column)
    	{
    		printf("ERROE: Not a square matrix!\n");
    		return;
    	}
    
    	int i, j, k, m;
    	int size;
    	const int N = MatrixRow(A);
    	ComplexType temp;
    
    	Matrix a, b;
    	InitComplexMatrix(&a, N, 1);
    	InitComplexMatrix(&b, N, 1);
    	size = MatrixSize(A);
    	if (MatrixSize(Q) != size)
    	{
    		DestroyComplexMatrix(Q);                   
    		InitComplexMatrix(Q, A->row, A->column); 
    	}
    
    	if (MatrixSize(R) != size)
    	{
    		DestroyComplexMatrix(R);
    		InitComplexMatrix(R, A->row, A->column);
    	}
    
    	for (j = 0; j < N; ++j)
    	{
    		for (i = 0; i < N; ++i)
    		{
    			a.arrayComplex[i] = b.arrayComplex[i] = A->arrayComplex[i * A->column + j];
    		}
    
    		for (k = 0; k < j; ++k)
    		{
    			R->arrayComplex[k * R->column + j]._Val[0] = 0;
    			R->arrayComplex[k * R->column + j]._Val[1] = 0;
    			for (m = 0; m < N; ++m)
    			{
    				R->arrayComplex[k * R->column + j] = AddComplex(R->arrayComplex[k * R->column + j], \
    					_Cmulcc(a.arrayComplex[m], Q->arrayComplex[m * Q->column + k]));
    			}
    
    			for (m = 0; m < N; ++m)
    			{
    				b.arrayComplex[m] = SubComplex(b.arrayComplex[m], _Cmulcc(R->arrayComplex[k * R->column + j], \
    					Q->arrayComplex[m * Q->column + k]));
    			}
    		}
    
    		temp = MatrixNorm2(&b);
    		R->arrayComplex[j * R->column + j] = temp;
    
    		for (i = 0; i < N; ++i)
    		{
    			Q->arrayComplex[i * Q->column + j] = DivComplex(b.arrayComplex[i], temp);
    		}
    	}
    	Matrix ComplexMatrixArray[] = { a, b };
    	int numDoubleMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
    	DestroyComplexMatrixArray(ComplexMatrixArray, numDoubleMatrixArray);
    }
    
    /*
      Complex Matrix Multiple: matrixC = matrixA * matrixB
    */
    void MatrixMulMatrix(const Matrix* matrixA, const Matrix* matrixB, Matrix* matrixC)
    {
    	if (IsNullComplexMatrix(matrixA) || IsNullComplexMatrix(matrixB))
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return;
    	}
    
    	else if (matrixA->column != matrixB->row)
    	{
    		puts("ERROE: An incompatable matrix!\n");
    		return;
    	}
    	else
    	{
    		int row_i, column_j, ij;
    		int indexA, indexB, indexC;
    		ComplexType tempCell_C;
    
    		for (row_i = 0; row_i < matrixC->row; ++row_i)
    		{
    			for (column_j = 0; column_j < matrixC->column; ++column_j)
    			{
    				tempCell_C._Val[0] = 0;
    				tempCell_C._Val[1] = 0;
    				for (ij = 0; ij < matrixA->column; ++ij)
    				{
    					indexA = row_i * matrixA->column + ij;
    					indexB = ij * matrixB->column + column_j;
    					tempCell_C = AddComplex(tempCell_C,
    						_Cmulcc(matrixA->arrayComplex[indexA], matrixB->arrayComplex[indexB]));
    				}
    				indexC = row_i * matrixC->column + column_j;
    				matrixC->arrayComplex[indexC] = tempCell_C;
    			}
    		}
    	}
    }
    
    /*
       eigen values
    */
    
    void EigenValue(const Matrix* matrix, Matrix* eigenvalue)
    {
    	const int NUM = 100; 
    
    	Matrix Q, R, temp;
    	InitComplexMatrix(&Q, matrix->row, matrix->column);
    	InitComplexMatrix(&R, matrix->row, matrix->column);
    	InitComplexMatrix(&temp, matrix->row, matrix->column);
    
    	if (IsNullComplexMatrix(matrix))
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return;
    	}
    
    	CopyMatrix(matrix, &temp);
    
    	for (int k = 0; k < NUM; ++k)
    	{
    		QR(&temp, &Q, &R);
    		MatrixMulMatrix(&R, &Q, &temp);
    	}
    	for (int k = 0; k < temp.row; ++k)
    	{
    		eigenvalue->arrayComplex[k] = temp.arrayComplex[k * temp.column + k];
    	}
    	Matrix ComplexMatrixArray[] = { R, Q, temp };
    	int numComplexMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
    	DestroyComplexMatrixArray(ComplexMatrixArray, numComplexMatrixArray);
    }
    
    /*
       eigen vectors
    */
    void EigenVector(const Matrix* matrix, const Matrix* eigenvalue, Matrix* eigenvector)
    {
    	if (IsNullComplexMatrix(matrix) || IsNullComplexMatrix(eigenvalue))
    	{
    		puts("ERROE: An invalid matrix!\n");
    		return;
    	}
    	int i, j, q;
    	int m;
    	int count;
    	int num = MatrixRow(matrix); 
    	ComplexType evalue;
    
    	ComplexType sum, midsum, mid;
    	Matrix temp;
    	InitComplexMatrix(&temp, matrix->row, matrix->column);
    
    	for (count = 0; count < num; ++count)
    	{
    		evalue = eigenvalue->arrayComplex[count];
    		CopyMatrix(matrix, &temp);
    		for (i = 0; i < temp.column; ++i)
    		{
    			temp.arrayComplex[i * temp.column + i] = SubComplex(temp.arrayComplex[i * temp.column + i], evalue);
    		}
    		for (i = 0; i < temp.row - 1; ++i)
    		{
    			mid._Val[0] = creal(temp.arrayComplex[i * temp.column + i]); 
    			mid._Val[1] = cimag(temp.arrayComplex[i * temp.column + i]);
    			for (j = i; j < temp.column; ++j)
    			{
    				temp.arrayComplex[i * temp.column + j] = DivComplex(temp.arrayComplex[i * temp.column + j], mid);
    			}
    
    			for (j = i + 1; j < temp.row; ++j)
    			{
    				mid = temp.arrayComplex[j * temp.column + i];
    				for (q = i; q < temp.column; ++q)
    				{
    					temp.arrayComplex[j * temp.column + q] = SubComplex(temp.arrayComplex[j * temp.column + q], \
    						_Cmulcc(temp.arrayComplex[i * temp.column + q], mid));
    
    				}
    			}
    		}
    		midsum._Val[0] = 1;
    		midsum._Val[1] = 0;
    		eigenvector->arrayComplex[(eigenvector->row - 1) * eigenvector->column + count]._Val[0] = 1;
    		eigenvector->arrayComplex[(eigenvector->row - 1) * eigenvector->column + count]._Val[1] = 0;
    		for (m = temp.row - 2; m >= 0; --m)
    		{
    			sum._Val[0] = 0; sum._Val[1] = 0;
    			for (j = m + 1; j < temp.column; ++j)
    			{
    				sum = AddComplex(sum, _Cmulcc(temp.arrayComplex[m * temp.column + j],
    					eigenvector->arrayComplex[j * eigenvector->column + count]));
    			}
    			sum = DivComplex(NegativeComplex(sum), temp.arrayComplex[m * temp.column + m]);
    			midsum = AddComplex(midsum, _Cmulcc(sum, sum));
    			eigenvector->arrayComplex[m * eigenvector->column + count] = sum;
    		}
    
    		midsum = csqrt(midsum);
    		for (i = 0; i < eigenvector->row; ++i)
    		{
    			eigenvector->arrayComplex[i * eigenvector->column + count] =
    				DivComplex(eigenvector->arrayComplex[i * eigenvector->column + count], midsum);
    		}
    	}
    	DestroyComplexMatrix(&temp);
    }
    
    int main(void)
    {
    	Matrix matrix;
    	InitComplexMatrix(&matrix, 3, 3);  // 3 * 3	
    	matrix.arrayComplex[0]._Val[0] = 2;  matrix.arrayComplex[0]._Val[1] = 2;
    	matrix.arrayComplex[1]._Val[0] = 3;  matrix.arrayComplex[1]._Val[1] = 7;
    	matrix.arrayComplex[2]._Val[0] = 3;  matrix.arrayComplex[2]._Val[1] = -6;
    	matrix.arrayComplex[3]._Val[0] = 6;  matrix.arrayComplex[3]._Val[1] = 2;
    	matrix.arrayComplex[4]._Val[0] = -7; matrix.arrayComplex[4]._Val[1] = 7;
    	matrix.arrayComplex[5]._Val[0] = 2;  matrix.arrayComplex[5]._Val[1] = 4;
    	matrix.arrayComplex[6]._Val[0] = -9;  matrix.arrayComplex[6]._Val[1] = 6;
    	matrix.arrayComplex[7]._Val[0] = 5;  matrix.arrayComplex[7]._Val[1] = -4;
    	matrix.arrayComplex[8]._Val[0] = -2;  matrix.arrayComplex[8]._Val[1] = -1;
    
    	Matrix EigenValues;
    	InitComplexMatrix(&EigenValues, MatrixRow(&matrix), 1);  // 3 * 1
    	EigenValue(&matrix, &EigenValues);
    
    	Matrix EigenVectors;
    	InitComplexMatrix(&EigenVectors, MatrixRow(&matrix), MatrixColumn(&matrix));  // 3 * 3
    	EigenVector(&matrix, &EigenValues, &EigenVectors);
    
    	for (int i = 0; i < EigenValues.row; i++)
    		for (int j = 0; j < EigenValues.column; j++)
    			printf("%lf, %lfi\n", creal(EigenValues.arrayComplex[i * EigenValues.column + j]), cimag(EigenValues.arrayComplex[i * EigenValues.column + j]));
    
    	printf("\n\n\n\n\n");
    
    	for (int index = 0; index < MatrixSize(&EigenVectors); index++)
    		printf("%lf, %lfi\n", creal(EigenVectors.arrayComplex[index]), cimag(EigenVectors.arrayComplex[index]));
    	// Free Array of Memories of Complex Matrice by Calling ComplexMatrixArray() 
    	Matrix ComplexMatrixArray[] = { matrix, EigenValues, EigenVectors };
    	int numComplexMatrixArray = sizeof(ComplexMatrixArray) / sizeof(Matrix);
    	DestroyComplexMatrixArray(ComplexMatrixArray, numComplexMatrixArray);
    	
    	return 0;
    

    测试结果我直接通过截图的形式给出,VS C(上)和matlab(下):
    在这里插入图片描述
    在这里插入图片描述
    PS
    (1) <2020.8.20补充> 上面看起来,二者的特征向量不一样,我刚开始也有点纳闷,后来经小伙伴@hamsterWalle的提示,发现matlab调用eig()函数求导的特征向量并未对向量进行归一化(归一化:特征向量/列向量各个元素的平方和为1),而在笔者的C程序中,对其进行了归一化,因此结果不一样,这里我做了两组实验:
    实验一:注意到前面给的求特征向量的函数EigenVector()中有两句程序

    midsum = AddComplex(midsum, _Cmulcc(sum, sum));
    
    midsum = csqrt(midsum);   // complex.h: csqrt()
    

    即是对特征向量进行归一化处理。此外,为了验证此函数求得的特征向量确实是归一化之后的,笔者在上面的demo中,main()函数之前加了一个归一化函数(注意是列向量,因此列向量的元素index与行向量有区别):

    /*
    	Normalization of Column Vector
    */
    void NormVector(const Matrix* EigenVector, Matrix* NormEigenVector)
    {
    	if (IsNullComplexMatrix(EigenVector))
    		return;
    	else if ((MatrixRow(EigenVector) != MatrixRow(NormEigenVector)) || (MatrixColumn(EigenVector) != MatrixColumn(NormEigenVector)))
    		return;
    	else
    	{
    		ComplexType exp = { 2, 0 };  // initiate a complex number
    		for (int column_j = 0; column_j < EigenVector->column; ++column_j)
    		{
    			ComplexType magnitude; InitComplex(&magnitude);
    			for (int row_i = 0; row_i < EigenVector->row; row_i++)
    			{
    				int index = row_i * EigenVector->column + column_j;
    				magnitude = AddComplex(magnitude, cpow(EigenVector->arrayComplex[index], exp));
    // OR: cpow() -----> _Cmulcc()
    //				magnitude = AddComplex(magnitude,  _Cmulcc(EigenVector->arrayComplex[index], EigenVector->arrayComplex[index]));
    			}
    
    			magnitude = csqrt(magnitude);
    
    			for (int row_i = 0; row_i < EigenVector->row; row_i++)
    			{
    				int index = row_i * NormEigenVector->column + column_j;
    				NormEigenVector->arrayComplex[index] = DivComplex(EigenVector->arrayComplex[index], magnitude);
    			}
    		}	
    	}
    }
    

    另外main函数内加上:

    	Matrix NormEigenVectors;
    	InitComplexMatrix(&NormEigenVectors, MatrixRow(&EigenVectors), MatrixColumn(&EigenVectors));  // 3 * 3
    	NormVector(&EigenVectors, &NormEigenVectors);
    // AND
    		printf("\n\n\n\n\n");
    	for (int index = 0; index < MatrixSize(&NormEigenVectors); index++)
    		printf("%lf, %lfi\n", creal(NormEigenVectors.arrayComplex[index]), cimag(NormEigenVectors.arrayComplex[index]));
    

    此外,如果要释放NormEIgenVectors矩阵的内存(其实是结构体内部指针的内存),还需要加一笔:

    Matrix ComplexMatrixArray[] = { matrix, EigenValues, EigenVectors, NormEigenVectors };
    

    最后编译结果完全一样,VS C结果如下,说明EigenVector()求得的特征向量已经归一化了,不需要单独写一个函数归一化:
    在这里插入图片描述
    实验二:为了验证上面VS和matlab求得的特征向量结果本质完全一样,只是matlab的eig()函数没有对特征向量归一化,笔者这里用matlab手动对上面matlab计算结果进行归一化,以特征值-12.060159-4.426611i对应的特征向量为例:
    在这里插入图片描述
    matlab结果截图如下:
    在这里插入图片描述
    在这里插入图片描述
    与前面VS求得的第一个特征值对应的特征向量(如下)完全一样:
    在这里插入图片描述
    (2) 关于matlab的eig()函数,若返回两个矢量,第一个是特征向量组成的矩阵,第二个是特征值组成的对角矩阵;若默认返回一个参数,则是一个特征值组成的列向量,此时matlab计算结果如下:
    在这里插入图片描述
    (3) 在求解特征向量函数中调用了前面章节没有涉及到的新的complex.h库函数—csqrt()函数和后面的归一化函数NormVector()中的cpow()函数,函数原型分别如下:

    _ACRTIMP _Dcomplex __cdecl csqrt(_In_ _Dcomplex _Z);
    
    _ACRTIMP _Dcomplex __cdecl cpow(_In_ _Dcomplex _X, _In_ _Dcomplex _Y);
    

    (4) 无论是matlab调用的eig()函数还是笔者编写的EigenValue()和EigenVector()函数,请读者注意维度,另外,求得的特征值和特征向量(矩阵)都是以列向量依次对应书写的,笔者都有标出。
    (5) <2020.8.20补充> 特征值若是实数,调用matlab的eig()函数会自动给特征值从小到大进行排序,读者若想实现跟matlab一样的效果,可以再写一个排序函数(快速/冒泡等排序算法),在调用完EigenValue()求得乱序的特征值之后,再调用该排序函数对特征值从小到大排序,然后再调用EigenValue()函数求特征向量,这样求得的特征值和特征向量会与matlab求得的顺序完全一致

    四、总结

    实数/复数矩阵的常见运算就到此为止,前后一起花了我大概四天的时间完成了这五篇连载博客的撰写,写demo、测试、总结等,希望能帮助到各位社区小伙伴!

    展开全文
  • 毕业多年,曾经有同事问我该如何理解特征值的意义?当时,实在羞愧,我一学数学的,真不知该如何回答。极力回想,也只能以“特征值的求法、步骤...bla...bla...”应付了事,答非所问,简直了得!这样的答案教科书里...

    1c5710317ced84c6f0bd76f0483edbec.png

    毕业多年,曾经有同事问我该如何理解特征值的意义?

    当时,实在羞愧,我一学数学的,真不知该如何回答。

    极力回想,也只能以“特征值的求法、步骤...bla...bla...”应付了事,

    答非所问,简直了得!

    这样的答案教科书里写得清清楚楚,网上Google/百度一大堆,

    人家问的是意义,如何理解其意义?

    直扣灵魂,

    我真的曾经理解过它的意义吗???

    招了吧,真没有!

    原在数学系时,教室里,对着黑板一堆密密麻麻的公式,我也是时常神游天外的主......

    考试前,为了避免挂科才熬夜突击,对着书本一一比划,至少要演算两到三张稿纸,才能勉强记住方法、步骤,哪还管得着它的意义?

    这种突击式的训练记忆,忘得也快,就像写代码一样,过一阵就忘了!

    课堂上,老师大多是照本宣科。

    当年,

    也许是自己知识阅历不够,很难理解其意义,

    也许是努力不够,被足球耽误了。

    也许是天赋所致,不能顿悟!

    ...

    总之,确定那时我肯定是没有理解它的意义的。

    不知道现在有多少学生还是一样?

    在学习一些抽象的数学工具时,代换三、四步之后就不知所云了,往往只能靠记忆强撑,而这种记忆最多维持一周,年轻时可能长点,后来,说忘就忘了......。

    有极少数天才,能在抽象世界里面一直转,抽啊抽,一直抽......并最终以此为业。

    而大多数人(99+%),一到毕业,就尴尬,因为真的不理解其意义,

    看似学了一些高深的数学知识,只会做题,不会运用,根本不理解公式指代符号的现实映射!进而职场上,其它方面训练缺失的短板逐渐显现后,囧是必然!

    我想,这不单是数学教育的问题,也是其它各方面可能会尴尬的本源:

    不理解意义

    好,扯远了,回到正题,来看灵魂之问:

    如何理解特征值的意义?

    最近才有些感悟,和大家分享一下。

    说到特征值,数学上,基本上是指矩阵的特征值

    说到矩阵,高等代数几乎一整本书都在讲它,最著名的数学软件叫Matlab,直译为矩阵实验室,足见其高深、复杂!

    而这么复杂混乱的东西确有一个特征值, 难道不奇怪?

    再说,矩阵到底有多复杂混乱?看数学公式体会一下吧:

    b106be4ab95f87ea26388918dae86f5f.png

    这是一堆数,每一个数字都可以在实数域内取值(正、负、零),或可以无限的延伸,联想到现在的大数据,还有什么东西不能由它表示呢?如果您相信万物皆数,这儿都可以说万物皆矩阵了,万物,能不复杂?

    另外,这一堆数既可以表示数据记录,还可以表示某种不知名的抽象运算(物理上叫算子),这样的数学运算,对某些对象集,确仅仅以固有的方式伸缩,且不管它是数据记录还是抽象运算,全都一样!

    如此混乱复杂! 确有本征!

    这不神奇吗?

    数学就是这样,抽象、高级、有理

    如果这样说感觉虚,那么先来看一下它精确(枯燥)的数学定义


    特征值

    设是一矩阵,是一维非零列向量,若存在一数,使得

    则称为的一个特征值,为的属于特征值的一个特征向量。

    展开

    , 即

    若把矩阵的每一行理解为一个基向量,则是表示基向量与该向量的内积() 等于。


    感觉公式真的很枯燥的话,就先跳过上面吧。

    下面我将从三个方面来试图阐释其意义,以便大家更好的理解。

    • 几何上
    • 医学上
    • 物理上

    (一)几何上

    如果把矩阵理解为一个坐标系(不一定直角的,也可能是严重变形的“尺子”),有一类向量叫特征向量,其中的任一个向量,在该矩阵上的投影,都只是自身固定的伸缩!

    如何理解投影呢?且拿三维来理解吧,一根射线在另外一个坐标系(矩阵)下的影子,其每一轴都会有投影分量,把所有分量组合还原成影子,跟其自身共线,影子射线的长度比值永远固定,这个比值就是特征值,简如下图。

    4c869b874bc0224d6080dd06a2647ce2.gif

    而该比值对这条直线上的所有向量都适应,即无论射线长短。那么有多少条这样的直线呢?维矩阵最多有条,每一条的比值(特征值)可能都不一样,有大有小,都代表这一维度的自身特征,故这里大、小意义就明显了。

    • 大可以理解为该维度上尺子的单位刻度大,比如表示一个单位刻度
    • 小可以理解为该维度上尺子的单位刻度小,比如表示一个单位刻度

    (二) 医学上

    如果把矩阵理解为中医祖传秘籍(乱不外传的),特征向量理解为秘方子(枸杞、百合、红花、童子尿...),特征值就是对该方子的用药量,温、热、寒不同方子特征值不一样, 这样也说得通,如下图!

    0a23fae744da3d7cf9fa2f1e4000a58f.png

    进一步,把西药制成品也类比为特征向量。比如新冠治疗中的瑞得西韦, 特征值就是该神药该服用多少?还有其它药方子,如莲花清瘟等,假设都能治疗新冠肺炎,但用量肯定是不一样的,即不同特征向量对应的特征值不一。

    27d378fa13c24ebececca50c8f8d53b0.png

    如此看来,特征值可理解为医学上药物用量的一个刻度,也是中西医互相密而不宣的沟通桥梁,正如下图的30aebc85888800b504ef3f02f4e8e247.png

    (三) 物理上

    “遇事不决,量子力学” 戏谑的表明了量子力学的高深、难懂!

    且看薛定谔方程的前半部分,就复杂得都让人头晕眼花.....

    9fddb12ec6233956290c535bf122b4cb.png

    物理学家把这种神操作统称为算子(因为给您解释不清楚~),是不是有点巫师作法、道士占卜的感觉?

    不同的是那帮巫师(物理学家),在圈内对不同公式符号都给出了互相认可的解释!

    例如:量子力学把世界看成是波动的,如果一个波函数经过一个量子变换后,它仍是一同一个波函数乘一个常量(如上图C)。

    ......

    再看矩阵,它不也就是一个算子吗?而且还是线性的,如此简单,so easy!

    大巫师(物理学家)牛!

    这样,特征值的意义又从矩阵的线性上升到非线性统一了。

    还是大巫师(物理学家)牛~

    总之,就是一段复杂的操作,统称为算子特征值也叫算子的本征值,台湾人习惯这样称呼,同一个意思,英文词源其实来自德语(自身的)。1a62d3b32f054a83b8aad56edcaab3ef.png

    本来很好理解的概念,几经"转手"之后就晦涩难懂了....

    遥想当年,若彼时能有这样的理解,就完美了!

    想记这点感悟也很久了,

    若有缘遇上,能给您带来一点点共鸣,便是满足506115fedf1569ee15161524bf4557d0.png

    最后附上特征值的求法,以便大家回忆。

    附:特征值求法


    特征多项式

    它是数域上的一个次多项式,若是复数域,必有个根。每一个根都是矩阵的一个特征值

    求特征值与特征向量方法步骤

    • 求出特征多项式的全部根
    • 把所求根代入方程组中求出一组基础解系,就得到属于相应特征值的线性无关的特征向量

    数学之水,更多免费干货、敬请关注详查~~

    c6356a5e5763e2ba756d576daad2fdc4.png

    423ed20820e598efe3e4a249d88ec541.png

    展开全文
  • 如何理解特征值复数的情况

    万次阅读 多人点赞 2019-10-03 13:49:58
    如何理解特征值复数的情况 特征值与特征向量 特征值定义为,若有Ax=λxAx=\lambda xAx=λx,则称xxx为AAA的特征向量,λ\lambdaλ为相应的特征值。这时我们可以发现,如果λλλ是实数,那么矩阵AAA对向量xxx的...

    如何理解特征值为复数的情况

    特征值与特征向量

    特征值可定义为,若有 A x = λ x Ax=\lambda x Ax=λx,则称 x x x A A A的特征向量, λ \lambda λ为相应的特征值。这时我们可以发现,如果 λ λ λ是实数,那么矩阵 A A A对向量 x x x的作用,恰恰是将向量 x x x放大或缩小了一定倍数,而方向不变。如果特征值绝对值大于1,那么 x x x的长度将被放大(方向可能改变),如果特征值绝对值小于1,那么 x x x的长度被缩小(方向也可能改变)

    特征向量在控制领域的意义

    由于答主本科阶段学的是自动化专业,接下来用控制领域的知识进一步阐述特征向量相关知识。特征向量是什么呢?就是系统状态的某个特殊值。而根据线性定理,总可以将系统状态拆成不同特征向量的线性组合,所以也可以认为特征向量是系统状态的一个分量。在(离散)控制系统中,若 ϕ \phi ϕ是系统离散域的状态转移矩阵,假如某个时刻系统的状态 x x x,恰好等于特征向量的倍数,那么每过一拍后,新的状态 x ( n + 1 ) = ϕ x ( n ) = λ x ( n ) x(n+1)=\phi x(n)=\lambda x(n) x(n+1)=ϕx(n)=λx(n)。如果 λ \lambda λ的绝对值大于1,那么 x x x中的值将随时间推移而不断变大,如果 λ \lambda λ绝对值小于1,那么 x x x中的值将随时间推移而不断衰减。系统的输出也包括在系统状态中,因此如果 λ \lambda λ绝对值大于1,系统就是发散的。这也正是Z变换中的结论。

    特征值为复数的情况

    题主问的是特征值为复数的情况。当 λ \lambda λ为复数时,从数学的角度看,可以认为状态 x x x每经过一个控制节拍(经过一拍意味着系统经过了一次状态转移,下同),不仅长度发生改变,还产生了一个复平面内的旋转。那么复平面上的旋转,与实际的物理意义有什么对应关系呢?在我看来,可以这么理解:若选取合适坐标系,复平面中的值,其实部表征了当前物理量的大小,虚部象征了物理量转化(如能量转移)的趋势。以理想的LC串联电路为例:

    电路工作过程为:若电感的感应电动势大于电容,则它给电容充能,否则电容给电感充能。由于没有能量的损耗,这个过程会持续到天荒地老。因此,若选取电容电压和电感电流构成一组系统状态,其离散系统状态转移矩阵的特征值就是“绝对值为1,而角度不为0”的一个复数。

    举个例子,若电容为 1 F 1F 1F,电感为 1 H 1H 1H系统的初始状态(电容电压,电感电流)为 ( 1 V 1V 1V 0 A 0A 0A),但在虚平面上应表示成( 1 V 1V 1V i A iA iA),其实部为( 1 V 1V 1V 0 A 0A 0A),后续将进一步揭示这种表现形式的原因。
    若选择合适的控制周期,使对应的离散系统状态转移矩阵 ϕ \phi ϕ的特征值为 i i i,即 e π 2 e^{\frac{\pi}{2}} e2π,那么经过一拍的时间后,系统状态变成了( i V iV iV − 1 A -1A 1A),取其实数部分,得物理状态( 0 V 0V 0V − 1 A -1A 1A),说明电容的能量全部转移到电感上了。再过一拍,系统状态为( − 1 V -1V 1V, − i A -iA iA), 取实部得物理状态( − 1 V -1V 1V, 0 A 0A 0A)。这样随着时间推移,接下来的系统状态为 ( − i V -iV iV, 1 A 1A 1A),( 1 V 1V 1V, i A iA iA) 等等,其实部为实际物理量,如此往复循环。

    如果我们把控制周期缩短,则特征值将会改变。若特征值为 e π 4 e^{\frac{\pi}{4}} e4π,那么我们就可以更清楚地看到系统状态变化的过程如下表所示。可以看出,电容电压、电感电流都如正弦振荡般往复循环。

    时间复平面上的状态值实数部分(电容电压、电感电流)角度表示
    0( 1 V 1V 1V i A iA iA)( 1 V 1V 1V 0 A 0A 0A)( e 0 e^{0} e0 e π 2 e^{\frac{\pi}{2}} e2π)
    1( 2 + 2 i 2 V \frac{\sqrt{2} +\sqrt{2}i }{2}V 22 +2 iV − 2 − 2 i 2 A -\frac{\sqrt{2} -\sqrt{2}i }{2}A 22 2 iA)( 2 2 V \frac{\sqrt{2}}{2}V 22 V − 2 2 A -\frac{\sqrt{2}}{2}A 22 A)( e π 4 e^{\frac{\pi}{4}} e4π e 3 π 4 e^{\frac{3\pi}{4}} e43π)
    2( i V iV iV − 1 A -1A 1A)( 0 V 0V 0V − 1 A -1A 1A)( e π 2 e^{\frac{\pi}{2}} e2π e π e^{\pi} eπ)
    3( − 2 − 2 i 2 V -\frac{\sqrt{2} -\sqrt{2}i }{2}V 22 2 iV − 2 + 2 i 2 A -\frac{\sqrt{2}+\sqrt{2}i }{2}A 22 +2 iA)( − 2 2 V -\frac{\sqrt{2}}{2}V 22 V − 2 2 A -\frac{\sqrt{2}}{2}A 22 A)( e 3 π 4 e^{\frac{3\pi}{4}} e43π e 5 π 4 e^{\frac{5\pi}{4}} e45π)
    4( − 1 V -1V 1V − i A -iA iA)( − 1 V -1V 1V 0 A 0A 0A)( e π e^{\pi} eπ e 3 π 2 e^{\frac{3\pi}{2}} e23π)
    5( − 2 + 2 i 2 V -\frac{\sqrt{2} +\sqrt{2}i }{2}V 22 +2 iV 2 − 2 i 2 A \frac{\sqrt{2} -\sqrt{2}i }{2}A 22 2 iA)( − 2 2 V -\frac{\sqrt{2}}{2}V 22 V 2 2 A ) \frac{\sqrt{2}}{2}A) 22 A)( e 5 π 4 e^{\frac{5\pi}{4}} e45π e 7 π 4 e^{\frac{7\pi}{4}} e47π)
    6( − i V -iV iV 1 A 1A 1A) 0 V 0V 0V 1 A 1A 1A)( e 3 π 2 e^{\frac{3\pi}{2}} e23π e 0 e^{0} e0)
    7( 2 − 2 i 2 V \frac{\sqrt{2} -\sqrt{2}i }{2}V 22 2 iV 2 + 2 i 2 A \frac{\sqrt{2}+\sqrt{2}i }{2}A 22 +2 iA)( 2 2 V \frac{\sqrt{2}}{2}V 22 V 2 2 A \frac{\sqrt{2}}{2}A 22 A)( e 7 π 4 e^{\frac{7\pi}{4}} e47π e π 4 e^{\frac{\pi}{4}} e4π)
    8( 1 V 1V 1V i A iA iA)( 1 V 1V 1V 0 A 0A 0A)( e 0 e^{0} e0 e π 2 e^{\frac{\pi}{2}} e2π)

    无论在复频域中或者是Z域中,但凡特征值为复数,都意味着系统有一对共轭极点。共轭极点的位置表示系统的固有频率及其幅度变化情况。为方便理解,上面以Z平面中的特点举例(采样时间T不同,Z变换的结果也不同,那么特征值也不一样)。若特征值都在单位圆上,因此系统输出会等幅振荡,这是上述的情况。如果特征值在单位圆外,振荡幅值将不断增加(比如自激振荡之类),特征值在单位圆内,振荡幅值会不断衰减(把上述电路串一个电阻)。

    上面讨论的是系统状态与特征向量方向相同的情况。但因为系统状态可以拆分成多个特征向量的线性组合,因而该论述对于任何的系统初始状态,都是适用的。

    最后,给大家推荐一本《线性代数的几何意义》,是它带给了我灵感。

    展开全文
  • 线性代数之——复数矩阵

    千次阅读 2019-11-29 14:00:52
    即使矩阵是实的,特征值和特征向量也经常会是复数。 1. 虚数回顾 虚数由实部和虚部组成,虚数相加时实部和实部相加,虚部和虚部相加,虚数相乘时则利用 i2=−1i^2=-1i2=−1。 在虚平面,虚数 3+2i3+2i3+2i 是位于...
  • 如何理解矩阵特征值的意义?

    千次阅读 2020-05-22 09:33:11
    如何理解矩阵特征值的意义? 毕业多年,曾经有同事问我该如何理解特征值的意义? 当时,实在羞愧,我一学数学的,真不知该如何回答。 ...... 多年后有点感悟了,就当还愿吧,本文将从几何、医学、物理三个视角试图...
  • 矩阵特征值和特征向量的性质

    千次阅读 2022-03-22 17:00:08
    3、设,则复数范围内,A恰有N个特征值 4、对于每个特征值,都有无穷个特征向量 证: 所以为满足为特征值的一个特征向量,则任意乘以一个非零数k,则任然为满足为特征值的一个特征向量 所以可以得出,为特征值时...
  • 从机器学习、量子计算、物理到许多数学和工程的问题,都可以通过找到一个矩阵特征值和特征向量来解决。根据定义(标量λ、向量v是特征值、特征向量A):视觉上,Av与特征向量v位于同一直线上。这里有些例子。然而,Ax...
  • 矩阵特征值

    2019-10-07 22:58:18
    ...如何理解矩阵特征值? 想要理解特征值,首先要理解矩阵相似。什么是矩阵相似呢?从定义角度就是:存在可逆矩阵P...
  • 本章开始,笔者将详细讲解实数和复数矩阵的模(范数)、协方差等运算,并给出函数和测试demo的C代码,以及与matlab计算的结果; 并且,本章也出现了我在编写函数、后期与matlab联调跑数据时候踩到的一个坑,读者稍微...
  • 在数列、微分方程、矩阵三个不同的领域中都见到了“特征值”这个名字,难道仅仅是因为他们都刻画了问题的特征吗?当时我想,一定不是的。或许这些特征值是同一个特征值。因为本人水平较低,这个...
  • 矩阵特征值和特征向量

    千次阅读 2021-03-27 14:47:19
    矩阵特征值和特征向量特征值的理解特征值 、特征向量、方阵的关系几何意义物理意义特征分解特征值和特征向量的含义的应用补充一点 特征值的理解 其实刚上大学的时候上的线性代数课上,也只是简单讲解了特征值和...
  • 稀疏矩阵稀疏矩阵定义:含有大量0元素的矩阵一个稀疏矩阵包括m*3项元素,其中m是原数组中非零项的个数。其第一列是行下标,第二列是列下标,第三列是非零项的储存一个浮点数需要8字节,一个下标值需要4字节,则...
  • //获取矩阵特征值 2*1 MatrixXd D_temp,D;//注意这里定义的MatrixXd里没有c D_temp = evals.real();//获取特征值实数部分 D = D_temp.asDiagonal(); cout ; cout ; //D cout ; //D MatrixXf::Index evalsMax; B....
  • 矩阵特征值分解

    千次阅读 2020-07-06 11:22:43
    特征值分解 物理意义: 矩阵可以表示一种变换; 特征向量表示矩阵变换的方向; 特征值表示矩阵变换在对应特征向量方向上的变换速度; 特征值与特征向量 ...从而,特征值与特征向量的定义式就是这样的:
  • 如果矩阵对某一个向量或某些向量只发生伸缩变换,不对这些向量产生旋转的效果,那么这些向量就称为这个矩阵的特征向量,伸缩的比例就是特征值。   实际上,上述的一段话既讲了矩阵变换特征值及特征向量的几何意义...
  • 矩阵论】特征值的估计(上下界和盖尔圆)

    千次阅读 多人点赞 2020-06-04 22:38:43
    因此,在实际工程计算上,面对高阶矩阵,我们常常通过计算特征值的范围来估计特征值,这个范围越小,所估计的特征值就越精确。以下,本文就特征值估计常用方法做以梳理,完整定理证明请参考西工大的《矩阵论》[1]。 ...
  • 如何理解矩阵特征值

    万次阅读 2016-05-19 12:31:15
    李浩 ,FPA蓝色 / EE。...特征值在很多领域应该都有自己的用途,它的物理意义到了本科高年级或者研究生阶段涉及到具体问题的时候就容易理解了,刚学线性代数的话,确实抽象。 ——————————————————
  • #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h>...//--------------------------这里是一些定义的...//定义矩阵元素的类型为matrixType typedef double matrixTy
  • 如何理解矩阵特征值

    千次阅读 2018-04-17 14:39:35
    对于运动而言,最重要的当然就是运动的速度和方向,那么(我后面会说明一下限制条件):特征值就是运动的速度特征向量就是运动的方向既然运动最重要的两方面都被描述了,特征值、特征向量自然可以称为运动(即矩阵)...
  • // 检验逆矩阵是否计算正确 write ( *,'(a)' ) '>> 检验结果是否为单位矩阵' x = matmul(matrix, x) do i = 1, n write(*,'(*(1x,f7.4,2x))') x(i,:) !// 检验结果是否为单位矩阵 end do deallocate( aa, ja, ia ) ...
  • 深入理解矩阵特征值和特征向量

    万次阅读 多人点赞 2019-09-16 16:29:40
    原 【数学基础】矩阵的特征向量、特征值及其含义 ...
  • 复习笔记的上篇在这里:辰晞:矩阵分析-期末复习笔记(上)​zhuanlan.zhihu.com目录:特征值,特征向量,相似 (Eigenvalues, eigenvectors, similarity)酉相似 & 酉等价 & 正规矩阵 (Unitary similarity &...
  • 我更愿意用:A^H=A来定义Hermite矩阵,这里A是任何一个n阶的复数矩阵,H表示共轭转置。你的这个问题有三个步骤,但是我不清楚你的意图是什么?你的问题(2)和(3)很容易回答,但是(1)我没有尝试怎么计算!!一定要用(1)...
  • 1、特征值是线性代数中的重要概念,设A是n阶方阵,如果存在数m和非零n维列向量x,使得Ax=mx成立,则称m是A的...当特征值复数时,就表示旋转90付。可用变换矩阵(0 -1 , 1 0)测试一下。 因为在复平面上,横轴表示实
  • c++使用Eigen库计算矩阵特征值

    千次阅读 2018-06-15 20:48:05
    下面的例子是计算矩阵特征值特征向量,并把特征向量矩阵实部和虚部分开。 #include &amp;lt;iostream&amp;gt; #include &amp;lt;Eigen/Dense&amp;gt; using namespace std; using namespace...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 5,681
精华内容 2,272
关键字:

复数矩阵的特征值定义

友情链接: trans.zip