C#，码海拾贝（20）——一般实矩阵的奇异值分解（Singular Value Decomposition）方法之C#源代码，《C#数值计算算法编程》源代码升级改进版
2023-09-11 14:15:48 时间
1 奇异值分解（Singular Value Decomposition）
奇异值分解（Singular Value Decomposition）是线性代数中一种重要的矩阵分解，奇异值分解则是特征分解在任意矩阵上的推广。在信号处理、统计学等领域有重要应用。
奇异值分解在某些方面与对称矩阵或Hermite矩阵基于特征向量的对角化类似。然而这两种矩阵分解尽管有其相关性，但还是有明显的不同。谱分析的基础是对称阵特征向量的分解，而奇异值分解则是谱分析理论在任意矩阵上的推广。
假设M是一个m×n阶矩阵，其中的元素全部属于域 K，也就是实数域或复数域。如此则存在一个分解使得
其中U是m×m阶酉矩阵；Σ是半正定m×n阶对角矩阵；而V*，即V的共轭转置，是n×n阶酉矩阵。这样的分解就称作M的奇异值分解。Σ对角线上的元素Σi，其中Σi即为M的奇异值。
常见的做法是为了奇异值由大而小排列。如此Σ便能由M唯一确定了。（虽然U和V仍然不能确定）
2 直观的解释
在矩阵M的奇异值分解中
·U的列(columns)组成一套对M的正交"输入"或"分析"的基向量。这些向量是MM*的特征向量。
·V的列(columns)组成一套对M的正交"输出"的基向量。这些向量是M*M的特征向量。
·Σ对角线上的元素是奇异值，可视为是在输入与输出间进行的标量的"膨胀控制"。这些是M*M及MM*的奇异值，并与U和V的列向量相对应。
奇异值和奇异向量, 以及他们与奇异值分解的关系
一个非负实数σ是M的一个奇异值仅当存在Km 的单位向量u和Kn的单位向量v如下 ：
其中向量u 和v分别为σ的左奇异向量和右奇异向量。
对于任意的奇异值分解，矩阵Σ的对角线上的元素等于M的奇异值。U和V的列分别是奇异值中的左、右奇异向量。因此，上述定理表明：
（1）一个m × n的矩阵至多有 p = min(m,n)个不同的奇异值；
（2）总是可以找到在Km 的一个正交基U，组成M的左奇异向量；
（3）总是可以找到和Kn的一个正交基V，组成M的右奇异向量。
如果一个奇异值中可以找到两个左（或右）奇异向量是线性相关的，则称为退化。
非退化的奇异值具有唯一的左、右奇异向量，取决于所乘的单位相位因子eiφ（根据实际信号）。因此，如果M的所有奇异值都是非退化且非零，则它的奇异值分解是唯一的，因为U中的一列要乘以一个单位相位因子且同时V中相应的列也要乘以同一个相位因子。
根据定义，退化的奇异值具有不唯一的奇异向量。因为，如果u1和u2为奇异值σ的两个左奇异向量，则两个向量的任意规范线性组合也是奇异值σ一个左奇异向量，类似的，右奇异向量也具有相同的性质。因此，如果M 具有退化的奇异值，则它的奇异值分解是不唯一的。
3 几何意义
因为U 和V 向量都是单位化的向量, 我们知道U的列向量u1,...,um组成了K空间的一组标准正交基。同样，V的列向量v1,...,vn也组成了K空间的一组标准正交基(根据向量空间的标准点积法则)。
线性变换T：即K → K，把向量Nx变换为Mx。考虑到这些标准正交基，这个变换描述起来就很简单了：T(vi) = σi ui, for i = 1,...,min(m,n)，其中σi 是对角阵Σ中的第i个元素；当i > min(m,n)时，T(vi) = 0。
这样，SVD理论的几何意义就可以做如下的归纳：对于每一个线性映射T: K → K，T把K的第i个基向量映射为K的第i个基向量的非负倍数,然后将余下的基向量映射为零向量。对照这些基向量，映射T就可以表示为一个非负对角阵。
4 应用
4.1 求伪逆
奇异值分解可以被用来计算矩阵的伪逆。若矩阵 M 的奇异值分解为 ，那么 M 的伪逆为。
其中 是 的伪逆，并将其主对角线上每个非零元素都求倒数之后再转置得到的。求伪逆通常可以用来求解线性最小平方、最小二乘法问题。
4.2 平行奇异值
把频率选择性衰落信道进行分解。
4.3 矩阵近似值
奇异值分解在统计中的主要应用为主成分分析（PCA），一种数据分析方法，用来找出大量数据中所隐含的“模式”，它可以用在模式识别，数据压缩等方面。PCA算法的作用是把数据集映射到低维空间中去。 数据集的特征值（在SVD中用奇异值表征）按照重要性排列，降维的过程就是舍弃不重要的特征向量的过程，而剩下的特征向量组成的空间即为降维后的空间。
5 一般实矩阵的奇异值分解 C# 源程序
5.1 主函数
using System;
namespace Zhou.CSharp.Algorithm
{
/// <summary>
/// 矩阵类
/// 作者：周长发
/// 改进：深度混淆
/// https://blog.csdn.net/beijinghorn
/// </summary>
public partial class Matrix
{
/// <summary>
/// 一般实矩阵的奇异值分解，分解成功后，原矩阵对角线元素就是矩阵的奇异值
/// </summary>
/// <param name="src">源矩阵</param>
/// <param name="mtxU">分解后的U矩阵</param>
/// <param name="mtxV">分解后的V矩阵</param>
/// <param name="eps">计算精度</param>
/// <returns>求解是否成功</returns>
public static bool SplitUV(Matrix src, Matrix mtxU, Matrix mtxV, double eps)
{
int i, j, k, u, it, uu, kk, ix, iy, mm, nn, iz, m1, ks;
double d, dd, t, sm, sm1, em1, sk, ek, b, c, shh;
double[] fg = new double[2];
double[] cs = new double[2];
int m = src.Rows;
int n = src.Columns;
// 初始化U, V矩阵
if (mtxU.Init(m, m) == false || mtxV.Init(n, n) == false)
{
return false;
}
// 临时缓冲区
int ka = Math.Max(m, n) + 1;
double[] s = new double[ka];
double[] e = new double[ka];
double[] w = new double[ka];
// 指定迭代次数为60
it = 60;
k = n;
if (m - 1 < n)
k = m - 1;
u = m;
if (n - 2 < m)
{
u = n - 2;
}
if (u < 0)
{
u = 0;
}
// 循环迭代计算
uu = k;
if (u > k)
{
uu = u;
}
if (uu >= 1)
{
for (kk = 1; kk <= uu; kk++)
{
if (kk <= k)
{
d = 0.0;
for (i = kk; i <= m; i++)
{
ix = (i - 1) * n + kk - 1;
d = d + src[ix] * src[ix];
}
s[kk - 1] = Math.Sqrt(d);
if (Math.Abs(s[kk - 1]) > float.Epsilon) // s[kk - 1] != 0.0)
{
ix = (kk - 1) * n + kk - 1;
if (Math.Abs(src[ix]) > float.Epsilon)//[ix] != 0.0)
{
s[kk - 1] = Math.Abs(s[kk - 1]);
if (src[ix] < 0.0)
{
s[kk - 1] = -s[kk - 1];
}
}
for (i = kk; i <= m; i++)
{
iy = (i - 1) * n + kk - 1;
src[iy] = src[iy] / s[kk - 1];
}
src[ix] = 1.0 + src[ix];
}
s[kk - 1] = -s[kk - 1];
}
if (n >= kk + 1)
{
for (j = kk + 1; j <= n; j++)
{
if ((kk <= k) && Math.Abs(s[kk - 1]) > float.Epsilon)//(s[kk - 1] != 0.0))
{
d = 0.0;
for (i = kk; i <= m; i++)
{
ix = (i - 1) * n + kk - 1;
iy = (i - 1) * n + j - 1;
d = d + src[ix] * src[iy];
}
d = -d / src[(kk - 1) * n + kk - 1];
for (i = kk; i <= m; i++)
{
ix = (i - 1) * n + j - 1;
iy = (i - 1) * n + kk - 1;
src[ix] = src[ix] + d * src[iy];
}
}
e[j - 1] = src[(kk - 1) * n + j - 1];
}
}
if (kk <= k)
{
for (i = kk; i <= m; i++)
{
ix = (i - 1) * m + kk - 1;
iy = (i - 1) * n + kk - 1;
mtxU[ix] = src[iy];
}
}
if (kk <= u)
{
d = 0.0;
for (i = kk + 1; i <= n; i++)
{
d = d + e[i - 1] * e[i - 1];
}
e[kk - 1] = Math.Sqrt(d);
if (Math.Abs(e[kk - 1]) > float.Epsilon)//e[kk - 1] != 0.0)
{
if (Math.Abs(e[kk]) > float.Epsilon)//e[kk] != 0.0)
{
e[kk - 1] = Math.Abs(e[kk - 1]);
if (e[kk] < 0.0)
{
e[kk - 1] = -e[kk - 1];
}
}
for (i = kk + 1; i <= n; i++)
{
e[i - 1] = e[i - 1] / e[kk - 1];
}
e[kk] = 1.0 + e[kk];
}
e[kk - 1] = -e[kk - 1];
if ((kk + 1 <= m) && Math.Abs(e[kk - 1]) > float.Epsilon)//(e[kk - 1] != 0.0))
{
for (i = kk + 1; i <= m; i++)
{
w[i - 1] = 0.0;
}
for (j = kk + 1; j <= n; j++)
{
for (i = kk + 1; i <= m; i++)
{
w[i - 1] = w[i - 1] + e[j - 1] * src[(i - 1) * n + j - 1];
}
}
for (j = kk + 1; j <= n; j++)
{
for (i = kk + 1; i <= m; i++)
{
ix = (i - 1) * n + j - 1;
src[ix] = src[ix] - w[i - 1] * e[j - 1] / e[kk];
}
}
}
for (i = kk + 1; i <= n; i++)
{
mtxV[(i - 1) * n + kk - 1] = e[i - 1];
}
}
}
}
mm = n;
if (m + 1 < n)
{
mm = m + 1;
}
if (k < n)
{
s[k] = src[k * n + k];
}
if (m < mm)
{
s[mm - 1] = 0.0;
}
if (u + 1 < mm)
{
e[u] = src[u * n + mm - 1];
}
e[mm - 1] = 0.0;
nn = m;
if (m > n)
{
nn = n;
}
if (nn >= k + 1)
{
for (j = k + 1; j <= nn; j++)
{
for (i = 1; i <= m; i++)
{
mtxU[(i - 1) * m + j - 1] = 0.0;
}
mtxU[(j - 1) * m + j - 1] = 1.0;
}
}
if (k >= 1)
{
for (uu = 1; uu <= k; uu++)
{
kk = k - uu + 1;
iz = (kk - 1) * m + kk - 1;
if (Math.Abs(s[kk - 1]) > float.Epsilon)//s[kk - 1] != 0.0)
{
if (nn >= kk + 1)
{
for (j = kk + 1; j <= nn; j++)
{
d = 0.0;
for (i = kk; i <= m; i++)
{
ix = (i - 1) * m + kk - 1;
iy = (i - 1) * m + j - 1;
d = d + mtxU[ix] * mtxU[iy] / mtxU[iz];
}
d = -d;
for (i = kk; i <= m; i++)
{
ix = (i - 1) * m + j - 1;
iy = (i - 1) * m + kk - 1;
mtxU[ix] = mtxU[ix] + d * mtxU[iy];
}
}
}
for (i = kk; i <= m; i++)
{
ix = (i - 1) * m + kk - 1;
mtxU[ix] = -mtxU[ix];
}
mtxU[iz] = 1.0 + mtxU[iz];
if (kk - 1 >= 1)
{
for (i = 1; i <= kk - 1; i++)
{
mtxU[(i - 1) * m + kk - 1] = 0.0;
}
}
}
else
{
for (i = 1; i <= m; i++)
{
mtxU[(i - 1) * m + kk - 1] = 0.0;
}
mtxU[(kk - 1) * m + kk - 1] = 1.0;
}
}
}
for (uu = 1; uu <= n; uu++)
{
kk = n - uu + 1;
iz = kk * n + kk - 1;
if ((kk <= u) && Math.Abs(e[kk - 1]) > float.Epsilon)
{
for (j = kk + 1; j <= n; j++)
{
d = 0.0;
for (i = kk + 1; i <= n; i++)
{
ix = (i - 1) * n + kk - 1;
iy = (i - 1) * n + j - 1;
d = d + mtxV[ix] * mtxV[iy] / mtxV[iz];
}
d = -d;
for (i = kk + 1; i <= n; i++)
{
ix = (i - 1) * n + j - 1;
iy = (i - 1) * n + kk - 1;
mtxV[ix] = mtxV[ix] + d * mtxV[iy];
}
}
}
for (i = 1; i <= n; i++)
{
mtxV[(i - 1) * n + kk - 1] = 0.0;
}
mtxV[iz - n] = 1.0;
}
for (i = 1; i <= m; i++)
{
for (j = 1; j <= n; j++)
{
src[(i - 1) * n + j - 1] = 0.0;
}
}
m1 = mm;
it = 60;
while (true)
{
if (mm == 0)
{
ppp(src.elements, e, s, mtxV.elements, m, n);
return true;
}
if (it == 0)
{
ppp(src.elements, e, s, mtxV.elements, m, n);
return false;
}
kk = mm - 1;
while ((kk != 0) && (Math.Abs(e[kk - 1]) > float.Epsilon)) //!= 0.0))
{
d = Math.Abs(s[kk - 1]) + Math.Abs(s[kk]);
dd = Math.Abs(e[kk - 1]);
if (dd > eps * d)
{
kk = kk - 1;
}
else
{
e[kk - 1] = 0.0;
}
}
if (kk == mm - 1)
{
kk = kk + 1;
if (s[kk - 1] < 0.0)
{
s[kk - 1] = -s[kk - 1];
for (i = 1; i <= n; i++)
{
ix = (i - 1) * n + kk - 1;
mtxV[ix] = -mtxV[ix];
}
}
while ((kk != m1) && (s[kk - 1] < s[kk]))
{
d = s[kk - 1];
s[kk - 1] = s[kk];
s[kk] = d;
if (kk < n)
{
for (i = 1; i <= n; i++)
{
ix = (i - 1) * n + kk - 1;
iy = (i - 1) * n + kk;
d = mtxV[ix];
mtxV[ix] = mtxV[iy];
mtxV[iy] = d;
}
}
if (kk < m)
{
for (i = 1; i <= m; i++)
{
ix = (i - 1) * m + kk - 1;
iy = (i - 1) * m + kk;
d = mtxU[ix];
mtxU[ix] = mtxU[iy];
mtxU[iy] = d;
}
}
kk = kk + 1;
}
it = 60;
mm = mm - 1;
}
else
{
ks = mm;
while ((ks > kk) && (Math.Abs(s[ks - 1]) > float.Epsilon))// != 0.0))
{
d = 0.0;
if (ks != mm)
{
d = d + Math.Abs(e[ks - 1]);
}
if (ks != kk + 1)
{
d = d + Math.Abs(e[ks - 2]);
}
dd = Math.Abs(s[ks - 1]);
if (dd > eps * d)
{
ks = ks - 1;
}
else
{
s[ks - 1] = 0.0;
}
}
if (ks == kk)
{
kk = kk + 1;
d = Math.Abs(s[mm - 1]);
t = Math.Abs(s[mm - 2]);
if (t > d)
{
d = t;
}
t = Math.Abs(e[mm - 2]);
if (t > d)
{
d = t;
}
t = Math.Abs(s[kk - 1]);
if (t > d)
{
d = t;
}
t = Math.Abs(e[kk - 1]);
if (t > d)
{
d = t;
}
sm = s[mm - 1] / d;
sm1 = s[mm - 2] / d;
em1 = e[mm - 2] / d;
sk = s[kk - 1] / d;
ek = e[kk - 1] / d;
b = ((sm1 + sm) * (sm1 - sm) + em1 * em1) / 2.0;
c = sm * em1;
c = c * c;
shh = 0.0;
//if ((b != 0.0) || (c != 0.0))
if (Math.Abs(b) > float.Epsilon || Math.Abs(c) > float.Epsilon)
{
shh = Math.Sqrt(b * b + c);
if (b < 0.0)
{
shh = -shh;
}
shh = c / (b + shh);
}
fg[0] = (sk + sm) * (sk - sm) - shh;
fg[1] = sk * ek;
for (i = kk; i <= mm - 1; i++)
{
sss(fg, cs);
if (i != kk)
{
e[i - 2] = fg[0];
}
fg[0] = cs[0] * s[i - 1] + cs[1] * e[i - 1];
e[i - 1] = cs[0] * e[i - 1] - cs[1] * s[i - 1];
fg[1] = cs[1] * s[i];
s[i] = cs[0] * s[i];
if ((Math.Abs(cs[0] - 1.0) > float.Epsilon) || Math.Abs(cs[1]) > float.Epsilon) //(cs[1] != 0.0))
{
for (j = 1; j <= n; j++)
{
ix = (j - 1) * n + i - 1;
iy = (j - 1) * n + i;
d = cs[0] * mtxV[ix] + cs[1] * mtxV[iy];
mtxV[iy] = -cs[1] * mtxV[ix] + cs[0] * mtxV[iy];
mtxV[ix] = d;
}
}
sss(fg, cs);
s[i - 1] = fg[0];
fg[0] = cs[0] * e[i - 1] + cs[1] * s[i];
s[i] = -cs[1] * e[i - 1] + cs[0] * s[i];
fg[1] = cs[1] * e[i];
e[i] = cs[0] * e[i];
if (i < m)
{
if ((cs[0] != 1.0) || Math.Abs(cs[1]) > float.Epsilon)//(cs[1] != 0.0))
{
for (j = 1; j <= m; j++)
{
ix = (j - 1) * m + i - 1;
iy = (j - 1) * m + i;
d = cs[0] * mtxU[ix] + cs[1] * mtxU[iy];
mtxU[iy] = -cs[1] * mtxU[ix] + cs[0] * mtxU[iy];
mtxU[ix] = d;
}
}
}
}
e[mm - 2] = fg[0];
it = it - 1;
}
else
{
if (ks == mm)
{
kk = kk + 1;
fg[1] = e[mm - 2];
e[mm - 2] = 0.0;
for (uu = kk; uu <= mm - 1; uu++)
{
i = mm + kk - uu - 1;
fg[0] = s[i - 1];
sss(fg, cs);
s[i - 1] = fg[0];
if (i != kk)
{
fg[1] = -cs[1] * e[i - 2];
e[i - 2] = cs[0] * e[i - 2];
}
if ((cs[0] != 1.0) || Math.Abs(cs[1]) > float.Epsilon)//(cs[1] != 0.0))
{
for (j = 1; j <= n; j++)
{
ix = (j - 1) * n + i - 1;
iy = (j - 1) * n + mm - 1;
d = cs[0] * mtxV[ix] + cs[1] * mtxV[iy];
mtxV[iy] = -cs[1] * mtxV[ix] + cs[0] * mtxV[iy];
mtxV[ix] = d;
}
}
}
}
else
{
kk = ks + 1;
fg[1] = e[kk - 2];
e[kk - 2] = 0.0;
for (i = kk; i <= mm; i++)
{
fg[0] = s[i - 1];
sss(fg, cs);
s[i - 1] = fg[0];
fg[1] = -cs[1] * e[i - 1];
e[i - 1] = cs[0] * e[i - 1];
if ((cs[0] != 1.0) || Math.Abs(cs[1]) > float.Epsilon)
{
for (j = 1; j <= m; j++)
{
ix = (j - 1) * m + i - 1;
iy = (j - 1) * m + kk - 2;
d = cs[0] * mtxU[ix] + cs[1] * mtxU[iy];
mtxU[iy] = -cs[1] * mtxU[ix] + cs[0] * mtxU[iy];
mtxU[ix] = d;
}
}
}
}
}
}
}
}
}
}
5.2 内部函数ppp，由SplitUV函数调用
POWER BY 315SOFT.COM
using System;
namespace Zhou.CSharp.Algorithm
{
/// <summary>
/// 矩阵类
/// 作者：周长发
/// 改进：深度混淆
/// https://blog.csdn.net/beijinghorn
/// </summary>
public partial class Matrix
{
/// <summary>
/// 内部函数，由SplitUV函数调用
/// </summary>
/// <param name="a"></param>
/// <param name="e"></param>
/// <param name="s"></param>
/// <param name="v"></param>
/// <param name="m"></param>
/// <param name="n"></param>
private static void ppp(double[] a, double[] e, double[] s, double[] v, int m, int n)
{
int i, j, p, q;
double d;
if (m >= n)
{
i = n;
}
else
{
i = m;
}
for (j = 1; j <= i - 1; j++)
{
a[(j - 1) * n + j - 1] = s[j - 1];
a[(j - 1) * n + j] = e[j - 1];
}
a[(i - 1) * n + i - 1] = s[i - 1];
if (m < n)
{
a[(i - 1) * n + i] = e[i - 1];
}
for (i = 1; i <= n - 1; i++)
{
for (j = i + 1; j <= n; j++)
{
p = (i - 1) * n + j - 1;
q = (j - 1) * n + i - 1;
d = v[p];
v[p] = v[q];
v[q] = d;
}
}
}
}
}
5.3 内部函数sss，由SplitUV函数调用
POWER BY TRUFFER.CN
using System;
namespace Zhou.CSharp.Algorithm
{
/// <summary>
/// 矩阵类
/// 作者：周长发
/// 改进：深度混淆
/// https://blog.csdn.net/beijinghorn
/// </summary>
public partial class Matrix
{
/// <summary>
/// 内部函数，由SplitUV函数调用
/// </summary>
/// <param name="fg"></param>
/// <param name="cs"></param>
private static void sss(double[] fg, double[] cs)
{
double r, d;
if (Math.Abs((Math.Abs(fg[0]) + Math.Abs(fg[1]))) < float.Epsilon)
{
cs[0] = 1.0;
cs[1] = 0.0;
d = 0.0;
}
else
{
d = Math.Sqrt(fg[0] * fg[0] + fg[1] * fg[1]);
if (Math.Abs(fg[0]) > Math.Abs(fg[1]))
{
d = Math.Abs(d);
if (fg[0] < 0.0)
{
d = -d;
}
}
if (Math.Abs(fg[1]) >= Math.Abs(fg[0]))
{
d = Math.Abs(d);
if (fg[1] < 0.0)
{
d = -d;
}
}
cs[0] = fg[0] / d;
cs[1] = fg[1] / d;
}
r = 1.0;
if (Math.Abs(fg[0]) > Math.Abs(fg[1]))
{
r = cs[1];
}
else if (Math.Abs(cs[0]) > float.Epsilon) //cs[0] != 0.0)
{
r = 1.0 / cs[0];
}
fg[0] = d;
fg[1] = r;
}
}
}
相关文章
- C# 系统应用之无标题窗体移动的两种方法
- C#图片压缩的实现方法
- c#通过socket判断服务器连接是否正常
- Win10系列：C#应用控件进阶9
- Camera Vision - video surveillance on C#
- 重学c#系列—— 反射深入一点点[三十三]
- C# 接口作用的深入理解
- C# DataTable.NewRow 方法
- C# 移除数组中重复数据
- C# DataTable.NewRow 方法
- C#小游戏之疯狂字母
- JS直接调用C#后台方法（ajax调用）
- C# 删除chart控件网格：两种方法
- C#里使用NPOI创建EXCEL文件的简单方法
- （60）C#里判断引用相等的方法
- C# 几种方法来复制的阵列
- .NET(C#) 读取Resource资源文件的方法
- C# prism 框架 MVVM框架 Prism系列之事件聚合器
- C#线程同步的几种方法
- C# 中使用using的三种方法