code: public function solve($B)
      {
          if ($B->getRowDimension() == $this->m) {
              if ($this->isFullRank()) {
                  // Copy right hand side
                  $nx = $B->getColumnDimension();
                  $X = $B->getArrayCopy();
                  // Compute Y = transpose(Q)*B
                  for ($k = 0; $k < $this->n; ++$k) {
                      for ($j = 0; $j < $nx; ++$j) {
                          $s = 0.0;
                          for ($i = $k; $i < $this->m; ++$i) {
                              $s += $this->QR[$i][$k] * $X[$i][$j];
                          }
                          $s = -$s / $this->QR[$k][$k];
                          for ($i = $k; $i < $this->m; ++$i) {
                              $X[$i][$j] += $s * $this->QR[$i][$k];
                          }
                      }
                  }
                  // Solve R*X = Y;
                  for ($k = $this->n - 1; $k >= 0; --$k) {
                      for ($j = 0; $j < $nx; ++$j) {
                          $X[$k][$j] /= $this->Rdiag[$k];
                      }
                      for ($i = 0; $i < $k; ++$i) {
                          for ($j = 0; $j < $nx; ++$j) {
                              $X[$i][$j] -= $X[$k][$j] * $this->QR[$i][$k];
                          }
                      }
                  }
                  $X = new Matrix($X);
  
                  return $X->getMatrix(0, $this->n - 1, 0, $nx);
              }
  
              throw new CalculationException(self::MATRIX_RANK_EXCEPTION);
          }
  
          throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
      }
docs: Least squares solution of A*X = B.
  
  @param Matrix $B a Matrix with as many rows as A and any number of columns
  
  @return Matrix matrix that minimizes the two norm of Q*R*X-B