В библиотеке OpenBLAS нет специальной функции для вычисления полярного разложения матрицы. Однако для этой цели может использоваться сингулярное разложение (SVD):

Формула полярного разложения: A = Q * P, где:

Q — ортогональная (или унитарная_ матрица

P — симметричная (или эрмитова) положительно определенная матрица.

Как можно рассчитать полярное разложение с помощью SVD.

Формула SVD: A = U * Σ * VT ==> (U * VT) * (V * Σ * VT), т.е. Q = U * VT и P = V * Σ * VT

Пример

//+------------------------------------------------------------------+
//| Тест функции PolarDecomposition: проверка Q ортогональна и Q*P=A |
//+------------------------------------------------------------------+
int TestPolarDecomposition(ulong size_m,ulong size_n)
  {
   matrix matrix_a=matrix::Random(size_m,size_n);
   matrix matrix_q,matrix_p;
   if(!PolarDecomposition(matrix_a,matrix_q,matrix_p))
      return(1000);
 
//--- проверка ортогональности
   matrix matrix_qqt=(size_m<=size_n) ? matrix_q@matrix_q.Transpose() : matrix_q.Transpose()@matrix_q;
   matrix matrix_i=matrix::Identity(matrix_qqt.Rows(),matrix_qqt.Cols());
   int    errors=(int)matrix_i.Compare(matrix_qqt,1e-11);
   printf("matrix Q is %sorthogonal",errors==0?"":"not ");
 
//--- тест полярного разложения Q * P = A
   matrix matrix_a2=matrix_q @ matrix_p;
   errors+=(int)matrix_a.Compare(matrix_a2,1e-11);
 
   printf("Polar Decomposition via SVD %s  matrix size %d x %d  errors=%d",errors==0?"passed":"failed",size_m,size_n,errors);
 
   return(errors);
  }
//+------------------------------------------------------------------+
//| Полярное разложение матрицы A = Q * P через сингулярное разлож.  |
//| Q — ортогональная, P — положит. определённая симметричная матрица|
//+------------------------------------------------------------------+
bool PolarDecomposition(matrixmatrix_a,matrixmatrix_q,matrixmatrix_p)
  {
//--- SVD
   matrix matrix_u,matrix_vt;
   vector vector_s;
   if(!matrix_a.SingularValueDecompositionDC(SVDZ_S,vector_s,matrix_u,matrix_vt))
      return(false);
 
//--- получаем диагональную матрицу
   ulong  diag=vector_s.Size();
   matrix matrix_s=matrix::Zeros(diag,diag);
   matrix_s.Diag(vector_s);
 
//--- Q = U * VT
   matrix_q=matrix_u @ matrix_vt;
//--- P = V * Σ * VT
   matrix_p=matrix_vt.Transpose() @ matrix_s @ matrix_vt;
 
   return(true);
  }