Calculator



合并了《Calculator V1.8.6.0（2013-10-31 10:33）》、《山寨版Matlab——Calculator V5.5beta(5.5版本更新 By2013/11/21 )2013-12-4 22:37》这几篇QQ日志。

如果A是一个m*n的矩阵，而B是一个p*q的矩阵，克罗内克积AⓧB则是一个mp*nq的分块矩阵。
矩阵的克罗内克积对应于线性映射的抽象张量积。并矢积是克罗内克积的特殊情况。克罗内克积是张量积的特殊形式，因此满足双线性与结合律，不符合交换律。
3行2列矩阵a与2行2列矩阵b的克罗内克积为6行4列的矩阵c
2行2列矩阵a与2行2列矩阵的克罗内克积为4行4列的矩阵c
如果A（m行n列），B（p行q列），C（n行c列），D（q行d列）是4个矩阵，且矩阵乘积AC（m行c列）和BD（p行d列）存在，那么：(AⓧB) (CⓧD)=

ACⓧBD=mp行cd列的矩阵=mp行nq列矩阵乘以nq行cd列矩阵
混合积的性质：它混合了通常的矩阵乘积和克罗内克积
AⓧB是可逆的当且仅当A（m行m列）和B（p行p列）是可逆的，其逆矩阵为：(AⓧB)^(-1)=A^(-1)ⓧB^(-1) =（mp行mp列）
如果A是n行n列矩阵，B是m行m列矩阵，I_k表示k行k列单位矩阵，那么我们可以定义克罗内克和⊕为：A⊕B=AⓧI_m+I_nⓧB=nm行nm列的矩阵
克罗内克积转置运算符合分配律：
如果A（m行n列），B（p行q列），(AⓧB)^(T)=A^(T)ⓧB^(T)=(mp行nq列)^(T)=nq行mp列的矩阵
rank(AⓧB)=rank(A)+rank(B)

功能:对矩阵进行Kronecker张量积.
格式:Kron(a,b)
说明:如果a是m*n的,b是p*q的,则kron(a,b)是大小(m*p)*(n*q)的矩阵
例子:
a=[1 2
3 1]
a =
[ 1.00000000000000    2.00000000000000
  3.00000000000000    1.00000000000000 ]
b=[0 3
2 1]
b =
[ 0.00000000000000    3.00000000000000
  2.00000000000000    1.00000000000000 ]
c=Kron(a,b)
c =
[ 0.00000000000000    3.00000000000000    0.00000000000000    6.00000000000000
  2.00000000000000    1.00000000000000    4.00000000000000    2.00000000000000
  0.00000000000000    9.00000000000000    0.00000000000000    3.00000000000000
  6.00000000000000    3.00000000000000    2.00000000000000    1.00000000000000 ]

a=[1
2
3]
a =
[ 1.00000000000000
  2.00000000000000
  3.00000000000000 ]
b=[4 5 6]
b =
[ 4.00000000000000    5.00000000000000    6.00000000000000 ]
c=Kron(a,b)
c =
[ 4.00000000000000    5.00000000000000    6.00000000000000
  8.00000000000000    10.0000000000000    12.0000000000000
  12.0000000000000    15.0000000000000    18.0000000000000 ]
d=mul(a,b)
d =
[ 4.00000000000000    5.00000000000000    6.00000000000000
  8.00000000000000    10.0000000000000    12.0000000000000
  12.0000000000000    15.0000000000000    18.0000000000000 ]

功能:求矩阵的2范数.谱范数.
格式:Cond2(a)
说明:a是任意维数的矩阵.本函数执行原理是,返回对称矩阵a'*a的最大特征值的平方根.其中a'是a的转置.
a=[1 -2
-3 4]
a =
[ 1.00000000000000   -2.00000000000000
 -3.00000000000000    4.00000000000000 ]
cond2(a)
ans =
[ 5.46498570421904 ]
矩阵A的2范数就是 A的转置乘以A矩阵特征根 最大值的开根号
如A={ 1 -2
-3 4 }
那么A的2范数就是msgbox (15+221^0.5)^0.5 '=5.46498570421904

范数是一种实泛函，满足非负性、齐次性、三角不等式
要理解矩阵的算子范数，首先要理解向量范数的内涵。矩阵的算子范数，是由向量范数导出的。
范数理论是矩阵分析的基础，度量向量之间的距离、求极限等都会用到范数，范数还在机器学习、模式识别领域有着广泛的应用。
矩阵范数反映了线性映射把一个向量映射为另一个向量，向量的"长度"缩放的比例。
由矩阵算子范数的定义形式可知，矩阵A把向量x映射成向量Ax，取其在向量x范数为1所构成的闭集下的向量Ax范数最大值作为矩阵A的范数，即矩阵对向量缩放的比例的上界，矩阵的算子范数是相容的。由几何意义可知，矩阵的算子范数必然大于等于矩阵谱半径（最大特征值的绝对值），矩阵算子范数对应一个取到向量Ax范数最大时的向量x方向，谱半径对应最大特征值下的特征向量的方向。而矩阵的奇异值分解SVD，分解成左右各一个酉阵，和拟对角矩阵，可以理解为对向量先作旋转、再缩放、最后再旋转，奇异值，就是缩放的比例，最大奇异值就是谱半径的推广，所以，矩阵算子范数大于等于矩阵的最大奇异值，酉阵在此算子范数的意义下，范数大于等于1。此外，不同的矩阵范数是等价的。

功能:方阵求逆,其采用的是LU分解算法
格式:inv(a)//a为方阵变量
a=[0.2368 0.2471 0.2568 1.2671
1.1161 0.1254 0.1397 0.149
0.1582 1.1675 0.1768 0.1871
0.1968 0.2071 1.2168 0.2271]
a =
[ 0.23680000000000     0.24710000000000     0.25680000000000     1.26710000000000
  1.11610000000000     0.12540000000000     0.13970000000000     0.14900000000000
  0.15820000000000     1.16750000000000     0.17680000000000     0.18710000000000
  0.19680000000000     0.20710000000000     1.21680000000000     0.22710000000000 ]
b=inv(a)
b =
[ -0.0859207504780      0.93794426823404     -0.0684372042645      -0.0796077151837
  -0.1055899132073      -0.0885243235004      0.90598255638825     -0.0991908105397
  -0.1270733117900      -0.1113511370480      -0.1169667064884      0.87842529094384
  0.85160581464323     -0.1354556628418      -0.1401825503018      -0.1438074804470 ]
ublas求矩阵的逆的输出及源代码：
mtxA矩阵[4,4]((0.2368,0.2471,0.2568,1.2671),(1.1161,0.1254,0.1397,0.149),(0.1582
,1.1675,0.1768,0.1871),(0.1968,0.2071,1.2168,0.2271))
mtxA逆矩阵[4,4]((-0.0859208,0.937944,-0.0684372,-0.0796077),(-0.10559,-0.0885243
,0.905983,-0.0991908),(-0.127073,-0.111351,-0.116967,0.878425),(0.851606,-0.1354
56,-0.140183,-0.143807))
mtxA逆矩阵的逆矩阵 [4,4]((0.2368,0.2471,0.2568,1.2671),(1.1161,0.1254,0.1397,0.1
49),(0.1582,1.1675,0.1768,0.1871),(0.1968,0.2071,1.2168,0.2271))
// REMEMBER to update "lu.hpp" header includes from boost-CVS
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/vector.hpp>

#include <iostream>

namespace ublas = boost::numeric::ublas;

//! \brief 求矩阵的逆
//! \param input[in] 被求矩阵
//! \param inverse[out] 逆矩阵
/* Matrix inversion routine.
Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */
template<class T>
bool InvertMatrix (const ublas::matrix<T>& input, ublas::matrix<T>& inverse) {
       using namespace boost::numeric::ublas;
       typedef permutation_matrix<std::size_t> pmatrix;
       // create a working copy of the input
       matrix<T> A(input);
       // create a permutation matrix for the LU-factorization
       pmatrix pm(A.size1());

       // perform LU-factorization
       int res = lu_factorize(A,pm);
       if( res != 0 ) return false;

       // create identity matrix of "inverse"
       inverse.assign(ublas::identity_matrix<T>(A.size1()));

       // backsubstitute to get the inverse
       lu_substitute(A, pm, inverse);

       return true;
}

int main()
{
       namespace ublas = boost::numeric::ublas;
       ublas::matrix<double> mTmp(4,4);
/*
    ' 原矩阵
    mtxA(1, 1) = 0.2368:    mtxA(1, 2) = 0.2471:    mtxA(1, 3) = 0.2568:    mtxA(1, 4) = 1.2671
    mtxA(2, 1) = 1.1161:    mtxA(2, 2) = 0.1254:    mtxA(2, 3) = 0.1397:    mtxA(2, 4) = 0.149
    mtxA(3, 1) = 0.1582:    mtxA(3, 2) = 1.1675:    mtxA(3, 3) = 0.1768:    mtxA(3, 4) = 0.1871
    mtxA(4, 1) = 0.1968:    mtxA(4, 2) = 0.2071:    mtxA(4, 3) = 1.2168:    mtxA(4, 4) = 0.2271
*/

#if 1
        //! \brief 求逆矩阵
        ublas::matrix<double> mtxA(4,4), MAInv(4,4);      
  mtxA(0, 0) = 0.2368;    mtxA(0, 1) = 0.2471;