//+------------------------------------------------------------------+
//|  Тест функции PseudoInverse: проверка A * A+ * A ≈ A             |
//+------------------------------------------------------------------+
int TestPseudoInverse(ulong size_m,ulong size_n)
  {
   matrix matrix_a=matrix::Random(size_m,size_n);
   matrix matrix_pi;
   if(!PseudoInverse(matrix_a,matrix_pi))
      return(1000);
 
//--- тест псевдо-инверсии A * A+ * A = A
   matrix matrix_a2=matrix_a @ matrix_pi @ matrix_a;
   int    errors=(int)matrix_a.Compare(matrix_a2,1e-11);
 
   printf("Pseudo Inverse via SVD %s  matrix size %d x %d  errors=%d",errors==0?"passed":"failed",size_m,size_n,errors);
 
   return(errors);
  }
//+------------------------------------------------------------------+
//| Вычисляет псевдообратную матрицу с помощью сингулярного разлож.  |
//| A+ = V * Σ⁻¹ * Uᵀ                                                |
//+------------------------------------------------------------------+
bool PseudoInverse(matrix& matrix_a,matrix& matrix_pi)
  {
//--- 1. SVD
   matrix matrix_u,matrix_vt;
   vector vector_s;
   if(!matrix_a.SingularValueDecompositionDC(SVDZ_S,vector_s,matrix_u,matrix_vt))
      return(false);
 
//--- 2. обратная диагональная матрица
   ulong  diag=vector_s.Size();
   matrix matrix_si=matrix::Zeros(diag,diag);
   for(ulong i=0; i<diag; i++)
     {
      if(vector_s[i]!=0.0)
         matrix_si[i][i]=1.0/vector_s[i];
     }
 
//--- 3. pseudo inverse  A+ = V * Σ+ * UT
   matrix_pi=matrix_vt.Transpose() @ matrix_si @ matrix_u.Transpose();
 
   return(true);
  }