Быстрое решение плотной линейной системы фиксированной размерности (N = 9), симметричной, положительно-полуопределенной

Как правило, лучше всего использовать существующую библиотеку, а не подход «по-своему», поскольку есть много утомительных деталей, на которые следует обратить внимание в стремлении к быстрой и стабильной численной реализации.

Вот несколько, чтобы вы начали:

Собственная библиотека (мои личные предпочтения):
http://eigen.tuxfamily.org/dox/QuickRefPage.html#QuickRef_Headers

Armadillo:
http://arma.sourceforge.net/

Ищите вокруг, и вы найдете много других.

Я бы порекомендовал разложение LU, особенно если «решить миллионы раз» действительно означает «решить один раз и применить к миллионам векторов». Вы создадите LU-декомпозицию, сохраните ее и примените прямую-обратную подстановку против стольких r.h.s. векторы, как вы хотите.

Это более устойчиво перед лицом округления, если вы используете поворот.

LU для симметричной полуопределенной матрицы не имеет особого смысла: вы уничтожаете красивое свойство ваших входных данных, выполняя ненужные операции.

Выбор между LLT или LDLT действительно зависит от числа условий ваших матриц и от того, как вы собираетесь обрабатывать крайние случаи. LDLT следует использовать только в том случае, если вы можете доказать статистически значимое улучшение точности или если надежность имеет первостепенное значение для вашего приложения.

(Без выборки ваших матриц трудно дать здравый совет, но я подозреваю, что при таком небольшом порядке N = 9 поворот маленьких диагональных членов к нижней части D действительно не нужен. Поэтому я бы начал с классической Cholesky и просто прервать разложение на множители, если термины diag станут слишком малыми по отношению к некоторому разумно выбранному допуску.

Cholesky довольно прост в коде, и если вы стремитесь к действительно быстрому коду, лучше реализовать его самостоятельно.

Как и другие выше, я рекомендую cholesky. Я обнаружил, что увеличение количества сложений, вычитаний и обращений к памяти означает, что LDLt медленнее, чем cholesky.

На самом деле существует несколько вариаций cholesky, и какая из них будет самой быстрой, зависит от того, какое представление вы выбрали для своих матриц. Я обычно использую представление в стиле фортрана, то есть матрица M является двойной * M, где M (i, j) равно m [i + dim * j]; для этого я считаю, что верхний треугольный холеский (немного) самый быстрый, то есть ищет верхний треугольный U с U ‘* U = M.

Для фиксированного, небольшого размера, безусловно, стоит подумать о написании версии, в которой нет циклов. Относительно простой способ сделать это — написать программу для этого. Насколько я помню, используя подпрограмму, которая работает с общим случаем в качестве шаблона, для написания программы, которая записывала бы конкретную версию с фиксированным размером, потребовалось только утро. Экономия может быть значительной. Например, моя общая версия занимает 0,47 секунды, чтобы выполнить миллион факторизаций 9×9, в то время как версия без петель занимает 0,17 секунды — эти временные интервалы работают на однопоточном компьютере с частотой 2,6 ГГц.

Чтобы показать, что это не главная задача, я включил источник такой программы ниже. Он включает в себя общую версию факторизации в качестве комментария. Я использовал этот код в обстоятельствах, когда матрицы не близки к единственному, и я считаю, что он хорошо работает там; однако это может быть слишком грубым для более деликатной работы.

/*  ----------------------------------------------------------------
**  to write fixed dimension ut cholesky routines
**  ----------------------------------------------------------------
*/
#include    <stdio.h>
#include    <stdlib.h>
#include    <math.h>
#include    <string.h>
#include    <strings.h>
/*  ----------------------------------------------------------------
*/
#if 0
static  inline  double  vec_dot_1_1( int dim, const double* x, const double* y)
{
double  d = 0.0;
while( --dim >= 0)
{   d += *x++ * *y++;
}
return d;
}

/*   ----------------------------------------------------------------
**  ut cholesky: solve U'*U = P for ut U in P (only ut of P accessed)
**  ----------------------------------------------------------------
*/

int mat_ut_cholesky( int dim, double* P)
{
int i, j;
double  d;
double* Ucoli;

for( Ucoli=P, i=0; i<dim; ++i, Ucoli+=dim)
{   /* U[i,i] = P[i,i] - Sum{ k<i | U[k,i]*U[k,i]} */
d = Ucoli[i] - vec_dot_1_1( i, Ucoli, Ucoli);
if ( d < 0.0)
{   return 0;
}
Ucoli[i] = sqrt( d);
d = 1.0/Ucoli[i];
for( j=i+1; j<dim; ++j)
{   /* U[i,j] = (P[i,j] - Sum{ k<i | U[k,i]*U[k,j]})/U[i,i] */
P[i+j*dim] = d*(P[i+j*dim] - vec_dot_1_1( i, Ucoli, P+j*dim));
}
}
return 1;
}
/*  ----------------------------------------------------------------
*/
#endif

/*  ----------------------------------------------------------------
**
**  ----------------------------------------------------------------
*/
static  void    write_ut_inner_step( int dim, int i, int off)
{
int j, k, l;
printf( "\td = 1.0/P[%d];\n", i+off);
for( j=i+1; j<dim; ++j)
{   k = i+j*dim;
printf( "\tP[%d] = d * ", k);
if ( i)
{   printf( "(P[%d]", k);
printf( " - (P[%d]*P[%d]", off, j*dim);
for( l=1; l<i; ++l)
{   printf( " + P[%d]*P[%d]", l+off, l+j*dim);
}
printf( "));");
}
else
{   printf( "P[%d];", k);
}
printf( "\n");
}
}

static  void    write_dot( int n, int off)
{
int i;
printf( "P[%d]*P[%d]", off, off);
for( i=1; i<n; ++i)
{   printf( "+P[%d]*P[%d]", off+i, off+i);
}
}

static  void    write_ut_outer_step( int dim, int i, int off)
{
printf( "\td = P[%d]", off+i);
if ( i)
{   printf( " - (");
write_dot( i, off);
printf( ")");
}
printf( ";\n");

printf( "\tif ( d <= 0.0)\n");
printf( "\t{\treturn 0;\n");
printf( "\t}\n");

printf( "\tP[%d] = sqrt( d);\n", i+off);
if ( i < dim-1)
{   write_ut_inner_step( dim, i, off);
}
}

static  void    write_ut_chol( int dim)
{
int i;
int off=0;
printf( "int\tut_chol_%.2d( double* P)\n", dim);
printf( "{\n");
printf( "double\td;\n");
for( i=0; i<dim; ++i)
{   write_ut_outer_step( dim, i, off);
printf( "\n");
off += dim;
}
printf( "\treturn 1;\n");
printf( "}\n");
}
/*  ----------------------------------------------------------------
*//*  ----------------------------------------------------------------
**
**  ----------------------------------------------------------------
*/
static  int read_args( int* dim, int argc, char** argv)
{
while( argc)
{   if ( strcmp( *argv, "-h") == 0)
{   return 0;
}
else if (strcmp( *argv, "-d") == 0)
{   --argc; ++argv;
*dim = atoi( (--argc, *argv++));
}
else
{   break;
}
}
return 1;
}

int main( int argc, char** argv)
{
int dim = 9;
if( read_args( &dim, --argc, ++argv))
{   write_ut_chol( dim);
}
else
{   fprintf( stderr, "usage: wchol (-d dim)? -- writes to stdout\n");
}
return EXIT_SUCCESS;
}
/*  ----------------------------------------------------------------
*/

Благодаря простоте использования, вы можете взять решатели Eigen только для сравнения. Для конкретного случая использования конкретный решатель может быть быстрее, хотя другой должен быть лучше. Для этого вы можете измерить время выполнения для каждого алгоритма только для выбора. После этого вы можете реализовать желаемый вариант (или найти существующий, который наилучшим образом соответствует вашим потребностям).