Replace Gram-Schmidt orthogonalisation with better algorithm

[mscroggs]

When elements are created, the polynomial basis is orthogonalised:

basix/cpp/basix/finite-element.cpp

Lines 41 to 62 in 0fc130f

	template <std::floating_point T>
	void orthogonalise(mdspan_t<T, 2>& wcoeffs)
	{
	for (std::size_t i = 0; i < wcoeffs.extent(0); ++i)
	{
	for (std::size_t j = 0; j < i; ++j)
	{
	T a = 0;
	for (std::size_t k = 0; k < wcoeffs.extent(1); ++k)
	a += wcoeffs(i, k) * wcoeffs(j, k);
	for (std::size_t k = 0; k < wcoeffs.extent(1); ++k)
	wcoeffs(i, k) -= a * wcoeffs(j, k);
	}
	T norm = 0;
	for (std::size_t k = 0; k < wcoeffs.extent(1); ++k)
	norm += wcoeffs(i, k) * wcoeffs(i, k);
	for (std::size_t k = 0; k < wcoeffs.extent(1); ++k)
	wcoeffs(i, k) /= std::sqrt(norm);
	}
	}

This orthogonalisation currently uses the Gram-Schmidt process. Instead, we should use a better algorithm (eg. Computing Q from the QR decomposition of A^T)

[ampdes]

The above function is now in math.h.

The Gram-Schmidt process constructs the Q matrix in-place whereas the R matrix can be constructed by collecting the a s (the a terms in above code snippet).

One alternative is the Householder reflection algorithm which constructs the R matrix in-place which is not sought here.

However, the current implementation is Classical Gram-Schmidt
https://github.com/FEniCS/basix/blob/44962386b6a5f06f64ec956cdd5737984df81f32/cpp/basix/math.h#L368C1-L391C1

can easily be changed to the numerically stable modified Gram-Schmidt:

  for (std::size_t i = start; i < wcoeffs.extent(0); ++i)
  {
    T norm = 0;
    for (std::size_t k = 0; k < wcoeffs.extent(1); ++k)
      norm += wcoeffs(i, k) * wcoeffs(i, k);
    if (std::abs(norm) < 4 * std::numeric_limits<T>::epsilon())
    {
      throw std::runtime_error(
          "Cannot orthogonalise the rows of a matrix with incomplete row rank");
    }

    for (std::size_t k = 0; k < wcoeffs.extent(1); ++k)
      wcoeffs(i, k) /= std::sqrt(norm);

    for (std::size_t j = i+1; j < wcoeffs.extent(0); ++j)
    {
      T a = 0;
      for (std::size_t k = 0; k < wcoeffs.extent(1); ++k)
        a += wcoeffs(i, k) * wcoeffs(j, k);
      for (std::size_t k = 0; k < wcoeffs.extent(1); ++k)
        wcoeffs(j, k) -= a * wcoeffs(i, k);
    }

  }

[mscroggs]

This orthogonalisation is now used in fewer places than it once was, so modified Gram-Schmidt should be good enough. Do you want to open a pull request to make the modification you suggested?

[ampdes]

This only affects the custom element creation?
For a PR, should this change be tested for test/test_custom_element.py?

[ampdes]

The existing implementation is neither the classical Gram Schmidt or the more commonly known modified Gram Schmidt.
A classical Gram Schmidt would use an auxiliary vector to store the orthogonalization result before scaling.
This version (existing implementation) appears in https://people.inf.ethz.ch/gander/papers/qrneu.pdf as a variant of modified Gram Schmidt.

None the less the commonly known modified Gram Schmidt is added by #785