Московский Государственный Технический Университет

имени Н.Э.Баумана

Кафедра САПР


Решение переопределенной системы линейных алгебраических уравнений методом SVD

Федорук В.Г.

Для решения системы линейных алгебраических уравнений (ЛАУ) вида AX=B, где A - прямоугольная матрица размерностью m строк на n столбцов, m>n, применение находит метод, основанный на SVD-разложении (Singular Value Decomposition) матрицы A:

A = U S VT

Здесь U - прямоугольная ортогональная по столбцам матрица размерностью m строк на n столбцов, V - квадратная ортогональная матрица размерностью n строк на n столбцов, S - диагональная матрица размерностью n строк на n столбцов, содержащая сингулярные значения матрицы A. Причем
UT U = VT V = E,

где E - единичная матрица размерностью n строк на n столбцов.

По найденному разложению матрицы A решение системы определяется следующим образом:

X = V diag{1/sii} UT B

Исходные тексты на языке СИ программ, реализующих SVD-разложение и решение на его базе систем ЛАУ, можно найти в литературе и в составе математических пакетов (например, http://www.netlib.org/clapack или http://alglib.sources.ru/matrixops/general).

Ниже приводится пример использования для решения системы ЛАУ функций из библиотеки GSL (GNU Scientific Library).

/*
* Решение системы ЛАУ A*X=B (кол-во уравнений > кол-во неизвестных) методом SVD
* (Singular Value Decomposition)
*
* Компиляция: gcc -o prg prg.c -L/usr/local/lib -lgsl -lgslcblas -lm
*/

 #include <stdio.h>
 #include <gsl/gsl_linalg.h>

int main (int argc, char *argv[]) {
/* Коэффициенты матрицы A[30x4] */
 double a_data[] = { 
 1.0, 0.1, 0.01, 0.001,
 1.0, 0.3, 0.09, 0.027,
 1.0, 0.5, 0.25, 0.125,
 1.0, 0.7, 0.49, 0.343,
 1.0, 0.9, 0.81, 0.729,
 1.0, 1.1, 1.21, 1.331,
 1.0, 1.3, 1.69, 2.197,
 1.0, 1.5, 2.25, 3.375,
 1.0, 1.7, 2.89, 4.913,
 1.0, 1.9, 3.61, 6.859,
 1.0, 2.1, 4.41, 9.261,
 1.0, 2.3, 5.29, 12.167,
 1.0, 2.5, 6.25, 15.625,
 1.0, 2.7, 7.29, 19.683,
 1.0, 2.9, 8.41, 24.389,
 1.0, 3.1, 9.61, 29.791,
 1.0, 3.3, 10.89, 35.937,
 1.0, 3.5, 12.25, 42.875,
 1.0, 3.7, 13.69, 50.653,
 1.0, 3.9, 15.21, 59.319,
 1.0, 4.1, 16.81, 68.921,
 1.0, 4.3, 18.49, 79.507,
 1.0, 4.5, 20.25, 91.125,
 1.0, 4.7, 22.09, 103.823,
 1.0, 4.9, 24.01, 117.649,
 1.0, 5.1, 26.01, 132.651,
 1.0, 5.3, 28.09, 148.877,
 1.0, 5.5, 30.25, 166.375,
 1.0, 5.7, 32.49, 185.193,
 1.0, 5.9, 34.81, 205.379
 };

 /* Элементы вектора правых частей B[30]*/
 double b_data[] = {
 0.0998334,
 0.29552,
 0.479426,
 0.644218,
 0.783327,
 0.891207,
 0.963558,
 0.997495,
 0.991665,
 0.9463,
 0.863209,
 0.745705,
 0.598472,
 0.42738,
 0.239249,
 0.0415807,
 -0.157746,
 -0.350783,
 -0.529836,
 -0.687766,
 -0.818277,
 -0.916166,
 -0.97753,
 -0.999923,
 -0.982453,
 -0.925815,
 -0.832267,
 -0.70554,
 -0.550686,
 -0.373877
 };

  /* Создание "представления" (view) для матрицы коэффициентов A из массива */
  gsl_matrix_view a = gsl_matrix_view_array (a_data, 30, 4);
  /* Теперь доступ к матрице через a.matrix */

  /* Создание "представления" для вектора B из массива */
  gsl_vector_view b = gsl_vector_view_array (b_data, 30);
  /* Теперь доступ к вектору через b.vector */

  gsl_matrix *v = gsl_matrix_alloc (4,4);	/* Вспомогательная матрица */
  gsl_vector *s = gsl_vector_alloc (4);	/* Вектор сингулярных значений */
  gsl_vector *work = gsl_vector_alloc (4);	/* Рабочий вектор */
  gsl_vector *x = gsl_vector_alloc (4);	/* Вектор X */

  /* Для доступа к элементам вектора можно использовать след. функции
  *	double gsl_vector_get (const gsl_vector * v, size_t i)
  *	void gsl_vector_set (gsl_vector * v, size_t i, double x)
  *
  *  Для доступа к элементам матрицы можно использовать след. функции
  *	double gsl_matrix_get (const gsl_matrix * m, size_t i, size_t j)
  *	void gsl_matrix_set (gsl_matrix * m, size_t i, size_t j, double x)
  *
  *  Индексация элементов - с 0
  */ 

  /* SVD - разложение матрицы A */
  gsl_linalg_SV_decomp (&a.matrix, v, s, work);
  /* Для коэффициентов A изменены значениями !!! */

  /* Решение системы ЛАУ */
  gsl_linalg_SV_solve (&a.matrix, v, s, &b.vector, x);

  /* Печать результатов */
  printf ("x = \n");
  gsl_vector_fprintf (stdout, x, "%g");

  gsl_vector_free (s);
  gsl_matrix_free (v);

    return 0;
  }