commit

Decomposition.cs

+using System;
+using System.Collections.Generic;
+using System.Text;
+
+namespace SVD
+{
+    static class Decomposition
+    {
+        public static bool outputSVD = false;
+        public static bool outputGolubKahan = false;
+
+        public static Matrix B22;
+        public static int p, q;
+
+        public static Tuple<double, double> Givens(double a, double b)
+        {
+            double sin, cos;
+            double tau;
+
+            if (b == 0)
+            {
+                cos = 1;
+                sin = 0;
+            }
+            else
+            {
+                if (Math.Abs(b) > Math.Abs(a))
+                {
+                    tau = -a / b;
+                    sin = 1 / Math.Sqrt(1 + tau * tau);
+                    cos = sin * tau;
+                }
+                else
+                {
+                    tau = -b / a;
+                    cos = 1 / Math.Sqrt(1 + tau * tau);
+                    sin = cos * tau;
+                }
+            }
+
+            return new Tuple<double, double>(cos, sin);
+        }
+
+        //Изменяет матрицу напрямую
+        public static Matrix RowRot(Matrix A, double cos, double sin)
+        {
+            int q = A.N;
+            double tau1, tau2;
+            
+            for (int i = 0; i < q; i++)
+            {
+                tau1 = A[0, i];
+                tau2 = A[1, i];
+
+                A[0, i] = cos * tau1 - sin * tau2;
+                A[1, i] = sin * tau1 + cos * tau2;
+            }
+
+            return A;
+        }
+
+        //Изменяет матрицу напрямую
+        public static Matrix ColRot(Matrix A, double cos, double sin)
+        {
+            int q = A.M;
+            double tau1, tau2;
+
+            for (int i = 0; i < q; i++)
+            {
+                tau1 = A[i, 0];
+                tau2 = A[i, 1];
+
+                A[i, 0] = cos * tau1 - sin * tau2;
+                A[i, 1] = sin * tau1 + cos * tau2;
+            }
+
+            return A;
+        }
+
+        //Изменяет матрицу B напрямую
+        public static Tuple<Matrix, Matrix, Matrix> GolubKahan(Matrix B)
+        {
+            //---------------------------------------//
+            //      Поиск собственного значения      //
+            //---------------------------------------//
+
+            //Обозначения по Вержбицкому
+            double p1, q2, pn1, qn1, pn, qn, tau;
+            double f, g;
+
+            int n = B.N;
+            
+            if (n > 2)
+            {
+                pn = B[n - 1, n - 1];
+                pn1 = B[n - 2, n - 2];
+
+                qn = B[n - 2, n - 1];
+                qn1 = B[n - 3, n - 2];
+
+                f = (qn1 * qn1 - qn * qn + pn1 * pn1 - pn * pn) / (2 * qn * pn1);
+
+                g = Math.Sqrt(1 + f * f);
+
+                tau = pn * pn + qn * qn + qn * pn1 * (f + g);
+            }
+            else
+            {
+                pn1 = B[0, 0];
+                pn = B[1, 1];
+                qn = B[0, 1];
+
+                tau = (1.0 / 2.0) * (pn1 * pn1 + pn * pn + qn * qn + Math.Sqrt(Math.Pow(pn1 * pn1 + pn * pn + qn * qn, 2) - 4 * pn1 * pn1 * pn * pn));
+            }
+
+
+            //---------------------------------------//
+
+            p1 = B[0, 0];
+            q2 = B[0, 1];
+
+            //Console.WriteLine("--------------");
+            //Console.WriteLine("Собственное значение: " + tau);
+            //Console.WriteLine("--------------");
+
+            double y = p1 * p1 - tau;
+            double z = p1 * q2;
+
+            //Console.WriteLine("Первый элемент в матрице BT * B: " + p1 * p1);
+            //Console.WriteLine("Второй элемент в матрице BT * B: " + p1 * q2);
+
+            Matrix U = Matrix.CreateIdentityMatrix(B.M, B.M);
+            Matrix V = Matrix.CreateIdentityMatrix(B.N, B.N);
+            Matrix U0;
+            Matrix V0;
+            double cos, sin;
+
+            Matrix matrix; //вспомогательная матрица
+
+            Matrix tempB = new Matrix(B); //удалить
+
+            //Console.WriteLine("=========== Начальное B22 =============");
+            //B.Output();
+
+            Tuple<double, double> cosSin;
+
+            for (int k = 0; k < n - 1; k++)
+            {
+                U0 = Matrix.CreateIdentityMatrix(B.M, B.M);
+                V0 = Matrix.CreateIdentityMatrix(B.N, B.N);
+
+                //--------------- Вычисление матрицы V0 -----------------//
+
+                cosSin = Givens(y, z);
+                cos = cosSin.Item1;
+                sin = cosSin.Item2;
+
+                matrix = V0.CreateSubMatrix(0, V0.M - 1, k, k + 1);
+                matrix = ColRot(matrix, cos, sin);
+                Matrix.AssignSubMatrix(V0, matrix, 0, k);
+                V = Matrix.Multiply(V, V0); //!!!
+                //Console.WriteLine("========== V0 ============");
+                //V0.Output();
+                //Console.WriteLine("========== V ============");
+                //V.Output();
+
+                matrix = B.CreateSubMatrix(0, B.M - 1, k, k + 1);
+                matrix = ColRot(matrix, cos, sin);
+                Matrix.AssignSubMatrix(B, matrix, 0, k);
+                //Console.WriteLine("========== B22 ============");
+                //B.Output();
+
+                //-------------------------------------------------------//
+
+                //--------------- Вычисление матрицы U0 -----------------//
+
+                y = B[k, k];
+                z = B[k + 1, k];
+
+                cosSin = Givens(y, z);
+                cos = cosSin.Item1;
+                sin = cosSin.Item2;
+
+                matrix = U0.CreateSubMatrix(k, k + 1, 0, U0.N - 1);
+                matrix = RowRot(matrix, cos, sin);
+                Matrix.AssignSubMatrix(U0, matrix, k, 0);
+                U = Matrix.Multiply(U0, U);
+                //Console.WriteLine("========== U0 ============");
+                //U0.Output();
+                //Console.WriteLine("========== U ============");
+                //U.Output();
+
+
+                matrix = B.CreateSubMatrix(k, k + 1, 0, B.N - 1);
+                matrix = RowRot(matrix, cos, sin);
+                Matrix.AssignSubMatrix(B, matrix, k, 0);
+                //Console.WriteLine("========== B22 ============");
+                //B.Output();
+
+                //-------------------------------------------------------//
+
+                if (k < n - 2)
+                {
+                    y = B[k, k + 1];
+                    z = B[k, k + 2];
+                }
+            }
+
+            //Console.WriteLine("############## Изменённая в Голубе-Кахане B22 ##############");
+            //B.Output();
+
+            //Console.WriteLine("================ Проверка ===============");
+            //Matrix.Multiply(Matrix.Multiply(U.CreateTransposeMatrix(), B), V.CreateTransposeMatrix()).Output();
+
+            //Console.WriteLine("======= Протестируем накопленные преобразования ======");
+            //Matrix.Multiply(Matrix.Multiply(U, tempB), V).Output();
+
+            return new Tuple<Matrix, Matrix, Matrix>(U, B, V);
+            //return new Tuple<Matrix, Matrix, Matrix>(U.CreateTransposeMatrix(), B, V.CreateTransposeMatrix());
+        }
+
+
+        public static void Decompose(Matrix B) //извлекает из матрицы B блок B22
+        {
+            //Console.WriteLine("======== Матрица на вход Decompose ========");
+            //B.Output();
+
+            int n = B.N;
+            int i = n - 2; //индекс предпоследнего диагонального элемента
+
+            int p; //размерность B11
+            int s, s1 = -1, s2 = -1; //вершины B22
+            bool firstPointIsFound = false;
+            int q; //размерность B33
+
+            while (i >= 0)
+            {
+                if (B[i, i + 1] != 0)
+                {
+                    if (!firstPointIsFound)
+                    {
+                        s2 = i + 1;
+                        firstPointIsFound = true;
+                    }
+                }
+                else
+                {
+                    if (firstPointIsFound)
+                    {
+                        s1 = i + 1;
+                        break;
+                    }
+                }
+
+                if (firstPointIsFound && i == 0)
+                {
+                    s1 = 0;
+                    break;
+                }
+
+                i--;
+            }
+
+            if (firstPointIsFound) //если матрица полностью диагональная
+            {
+                p = s1;
+                s = s2 - s1 + 1;
+                q = n - p - s;
+
+                Decomposition.B22 = B.CreateSubMatrix(p, p + s - 1, p, p + s - 1);
+            }
+            else
+            {
+                p = 0;
+                s = 0;
+                q = n;
+
+                Decomposition.B22 = B;
+            }
+
+            Decomposition.p = p;
+            Decomposition.q = q;
+        }
+
+        //Собирает матрицу B из блоков
+        public static Matrix Diag(Matrix A1, Matrix A2, Matrix A3)
+        {
+            int m = A1.M + A2.M + A3.M;
+
+            Matrix A = Matrix.CreateIdentityMatrix(m, m);
+
+            Matrix.AssignSubMatrix(A, A1, 0, 0);
+            Matrix.AssignSubMatrix(A, A2, A1.M, A1.M);
+            Matrix.AssignSubMatrix(A, A3, A1.M + A2.M, A1.M + A2.M);
+
+            return A;
+        }
+
+        public static Tuple<Matrix, Matrix, Matrix> SVD(Matrix A)
+        {
+            if (outputSVD)
+            {
+                Console.WriteLine("================================================================================================\n");
+                Console.WriteLine("                                   Сингулярное разложение                                       \n");
+                Console.WriteLine("================================================================================================\n");
+            }
+
+
+            double eps = 0.00000001;
+            int m = A.M;
+            int n = A.N;
+            int q = 0;
+            int p;
+            bool isZero = false; //некоторый диаг. элемент равен нулю
+
+            Matrix B = HouseHolder.CreateBidiagonalMatrix(A); //двухдиагонализация
+            Matrix U = HouseHolder.P.CreateTransposeMatrix();
+            Matrix V = HouseHolder.Q.CreateTransposeMatrix();
+
+            //Console.WriteLine("***** Проверка двухдиагонализации ******");
+
+            //U.Output();
+            //B.Output();
+            //V.Output();
+
+            //Console.WriteLine("*****************************************");
+
+            //Console.WriteLine("***** Проверка двухдиагонализации ******");
+            //B.Output();
+            //Matrix.Multiply(Matrix.Multiply(U, B), V).Output();
+
+
+            Matrix U0, V0; //для алгоритма 8.3.1
+            Matrix V1; //для обнуления наддиагональных элементов
+            Tuple<double, double> cosSin;
+
+            Matrix B22;
+
+            if (outputSVD)
+            {
+                Console.WriteLine("\n=========================  1. Двухдиагонализация  ===============================\n");
+                B.Output();
+            }
+
+
+
+            while (q != n)
+            {
+                //Зануляем наддиагональные элементы, значения которых меньше eps
+                for (int i = 0; i < n - 1; i++)
+                {
+                    if (Math.Abs(B[i, i + 1]) <= eps)
+                    {
+                        B[i, i + 1] = 0;
+                    }
+                }
+
+                Decompose(B); //результат Decompose сохраняется в стат. переменных
+
+                p = Decomposition.p;
+                q = Decomposition.q;
+                B22 = Decomposition.B22;
+
+
+                if (q < n)
+                {
+                    if (outputSVD)
+                    {
+                        Console.WriteLine("\n=========================  2.1. Выделение блока матрицы для метода Голуба-Канаха  ===================================\n");
+                        B22.Output();
+                    }
+
+                    //--------------------------------------------------------------------//
+                    //         Обнуление наддиагональных элементов той же строки,         //
+                    //                    если диагональный равен нулю                    //
+                    //--------------------------------------------------------------------//
+
+                    isZero = false;
+
+                    for (int i = 0; i < B22.M - 1; i++) //до предпоследнего диаг. элемента
+                    {
+                        if (B22[i, i] == 0)
+                        {
+                            V1 = Matrix.CreateIdentityMatrix(B22.N, B22.N);
+
+                            isZero = true;
+
+                            for (int j = i; j < B22.N; j++)
+                            {
+                                cosSin = Givens(B22[i, i], B22[i, i + 1]);
+
+                                V1 = RowRot(V1, cosSin.Item1, cosSin.Item2);
+
+
+                                V = Matrix.Multiply(V, V1); //накопление V
+
+                                B22 = RowRot(B22, cosSin.Item1, cosSin.Item2);
+                            }
+                        }
+                    }
+
+                    //------------------------------------------------------//
+
+                    if (!isZero)
+                    {
+                        Tuple<Matrix, Matrix, Matrix> UBV = GolubKahan(B22);
+                        U0 = UBV.Item1;
+                        B22 = UBV.Item2;
+                        V0 = UBV.Item3;
+
+                        Matrix.AssignSubMatrix(B, B22, p, p);
+
+                        if (outputSVD)
+                        {
+                            Console.WriteLine("\n======================  2.2 Блок матрицы после применения метода Голуба-Канаха  ================================\n");
+                            B22.Output();
+
+                            Console.WriteLine("\n=========================================  2.3 Полученная матрица  =============================================\n");
+                            B.Output();
+                        }
+
+                        //Console.WriteLine("======= Протестируем накопленные преобразования ======");
+                        //Matrix.Multiply(Matrix.Multiply(U0.CreateTransposeMatrix(), B22), V0.CreateTransposeMatrix()).Output();
+
+                        //Console.WriteLine("================== Матрица B до умножения ==================");
+                        //B.Output();
+
+                        Matrix Ip = Matrix.CreateIdentityMatrix(p, p);
+                        Matrix Iqmn = Matrix.CreateIdentityMatrix(q + m - n, q + m - n);
+                        Matrix D1 = Diag(Ip, U0, Iqmn);
+
+                        Matrix Iq = Matrix.CreateIdentityMatrix(q, q);
+                        Matrix D2 = Diag(Ip, V0, Iq);
+
+                        //Console.WriteLine("======================== D1 =========================");
+                        //D1.Output();
+                        //Console.WriteLine("======================== D2 ===============================");
+                        //D2.Output();
+
+                        U = Matrix.Multiply(D1, U); //накопление U
+                        V = Matrix.Multiply(V, D2); //накопление V
+
+
+                        //B = Matrix.Multiply(Matrix.Multiply(D1, B), D2);
+
+                        //Console.WriteLine("\n=============================== Матрица B после умножения ===============================\n");
+                        //B.Output();
+                    }
+                }
+            }
+
+            for (int i = 0; i < n; i++)
+            {
+                if (B[i, i] < 0)
+                {
+                    Matrix I = Matrix.CreateIdentityMatrix(n, n);
+                    I[i, i] = -1;
+                    B = Matrix.Multiply(B, I);
+                    V = Matrix.Multiply(V, I);
+                }
+            }
+
+            List<double> singValues = new List<double>();
+
+            for (int i = 0; i < B.N; i++)
+            {
+                singValues.Add(B[i, i]);
+            }
+
+            singValues.Sort();
+
+            Console.WriteLine("Сингулярные значения:");
+
+            for (int i = 0; i < singValues.Count; i++)
+            {
+                Console.WriteLine("{0}: {1}", i + 1, singValues[i]);
+            }
+
+            Console.WriteLine("================================================================================================\n");
+            Console.WriteLine("                              Завершение сингулярного разложения                                \n");
+            Console.WriteLine("================================================================================================\n");
+
+            return new Tuple<Matrix, Matrix, Matrix>(U, B, V);
+        }
+    }
+}
\ No newline at end of file

Matrix.cs

@@ -121,7 +121,7 @@ namespace SVD
 
                 for (int j = 0; j < n; j++)
                 {
-                    Console.Write("{0, 18}", Math.Round(data[i, j], 11));
+                    Console.Write("{0, 18}", Math.Round(data[i, j], 8));
                 }
                 Console.WriteLine("{0, 18}", "║");
             }
@@ -175,6 +175,11 @@ namespace SVD
                 throw new InvalidOperationException("Determinant can be calculated only for square matrix.");
             }
 
+            if (this.N == 1)
+            {
+                return this[0, 0];
+            }
+
             if (this.N == 2)
             {
                 return this[0, 0] * this[1, 1] - this[0, 1] * this[1, 0];
@@ -217,24 +222,6 @@ namespace SVD
             return result;
         }
 
-        public Matrix Transp()
-        {
-            int m = this.M;
-            int n = this.n;
-
-            Matrix result = new Matrix(n, m);
-
-            for (int i = 0; i < result.M; i++)
-            {
-                for (int j = 0; j < result.N; j++)
-                {
-                    result[i, j] = this[j, i];
-                }
-            }
-
-            return result;
-        }
-
         public Matrix CreateSubMatrix(int a1, int a2, int b1, int b2)
         {
             Matrix result = new Matrix(a2 - a1 + 1, b2 - b1 + 1);
@@ -305,7 +292,7 @@ namespace SVD
 
             if (n == 1)
             {
-                Matrix vectorT = this.Transp();
+                Matrix vectorT = this.CreateTransposeMatrix();
 
                 double element = Matrix.Multiply(vectorT, this)[0, 0];
 
@@ -405,5 +392,57 @@ namespace SVD
                 throw new Exception();
             }
         }
+
+        public static Matrix CreateDifferenceMatrix(Matrix A, Matrix B)
+        {
+            if (A.M == B.M && A.N == B.N)
+            {
+                int m = A.M;
+                int n = A.N;
+
+                Matrix result = new Matrix(m, n);
+
+                for (int i = 0; i < m; i++)
+                {
+                    for (int j = 0; j < n; j++)
+                    {
+                        double value = A[i, j] - B[i, j];
+                        result[i, j] = value;
+
+                    }
+                }
+
+                return result;
+            }
+            else
+            {
+                throw new Exception();
+            }
+        }
+
+        public static Matrix CreateIdentityMatrix(int m, int n)
+        {
+            Matrix matrix = new Matrix(m, n);
+
+            matrix.ProcessFunctionOverData((i, j) => matrix[i, j] = i == j ? 1 : 0);
+
+            return matrix;
+        }
+
+        public static void AssignSubMatrix(Matrix A, Matrix subA, int m, int n)
+        {
+            int m1 = m;
+            int m2 = m + subA.M - 1;
+            int n1 = n;
+            int n2 = n + subA.N - 1;
+
+            for (int i = m1; i <= m2; i++)
+            {
+                for (int j = n1; j <= n2; j++)
+                {
+                    A[i, j] = subA[i - m1, j - n1];
+                }
+            }
+        }
     }
 }

Program.cs

@@ -8,43 +8,335 @@ namespace SVD
     {
         static void Main(string[] args)
         {
-            //Matrix C = new Matrix(4, 3);
+            HouseHolder.output = true; //включение вывода двухдиагонализации
+            //Decomposition.outputGolubKahan = true;
+            Decomposition.outputSVD = false;
+
+            #region [Example1] Пример из учебника в Вержбицком 
+
+            //Matrix A = new Matrix(4, 3);
+
+            //A[0, 0] = 2;
+            //A[0, 1] = -1;
+            //A[0, 2] = -2;
+
+            //A[1, 0] = 0;
+            //A[1, 1] = 3;
+            //A[1, 2] = 1;
+
+            //A[2, 0] = -1;