Add LU solver documentation

[mgovers]

• Add LU solver documentation
  • Rationale
  • Implementation

[mgovers]

@TonyXiang8787 this is ready for a first review round. I will spend some time tomorrow on going through all the documentation from start to finish myself as well to check if the structure still makes sense, but some feedback would be nice as well at this stage.

[mgovers]

this is a potential follow-up improvement to the LU solver implementation

[TonyXiang8787]

I do not understand the logic here? What do you mean twice? The first time is just the solution of the original formula $Ax=b$ right? Because you initialize $r=b$. So you do only one actual iterative refinement.

[mgovers]

we always do one actual iterative refinement, but in the limit where $\epsilon_{backwarderror} \to \infty$, iterative refinement is never needed. If the threshold is set to infinity, it should break after the first time it goes through the loop, but in the current implementation, it will always do another solve-and-refine step

[TonyXiang8787]

I still do not understand your point here. So we now do the following initialization

1. berr <- inf.
2. r = b

The loop is (ignore the max limit)

while (berr > epsilon_converge)

1. solve A*dx = r
2. update berr, and update x
3. update r

So

1. In the first iteration, we are just solving Ax=b, which is the original equation. The berr will be very big. This should not be seen as iterative refinement.
2. In the second iteration, we do one iterative refinement. And we get the new berr.
3. If now the new berr is within the threshold, the loop terminates.

We always do ONE iterative refinement. But you mentioned in the doc that "iterative refinement will always be twice".

[mgovers]

This is another potential follow-up improvement to the LU solver implementation:

Because matrix inversion is equivalent to solving this set of equations on a matrix once, it actually may be faster to pre-calculate and then left multiply each block in the above by $\underline{m}^{-1}$ iff the row contains more than one element (which is pretty much always the case), and the difference in performance increases as the amount of off-diagonal elements increases.

In addition, it may make the code simpler. (TODO: check that that solution would be numerically stable)

[TonyXiang8787]

This was the very first prototype of the LU Solver. It quickly proved to be extremely numerically not statable. Calculating inverse explicitly you lose already half of the precision. Applying that inverse to other off-diagonal blocks makes the result much less accurate and leads to divergence of outer loop (NRPF).

Have a look at the original PR, especially this commit

08ee7fb

[TonyXiang8787]

You can explicitly document now in the docs to state applying inverse of diagonal blocks will result in numerical instability.

[TonyXiang8787]

Another thing I have tried is to use EIGEN build-in function to do the triangular solver instead of writing them manually. However, the benchmark was not so convincing. For some reason the EIGEN triangular solver is much slower than manually written stuff. Maybe this can be an investigation to understand why.

2cfe09b

[mgovers]

This was the very first prototype of the LU Solver. It quickly proved to be extremely numerically not statable. Calculating inverse explicitly you lose already half of the precision. Applying that inverse to other off-diagonal blocks makes the result much less accurate and leads to divergence of outer loop (NRPF).

Have a look at the original PR, especially this commit

08ee7fb

You can explicitly document now in the docs to state applying inverse of diagonal blocks will result in numerical instability.

Right. Thanks, that's useful input! Done.

Another thing I have tried is to use EIGEN build-in function to do the triangular solver instead of writing them manually. However, the benchmark was not so convincing. For some reason the EIGEN triangular solver is much slower than manually written stuff. Maybe this can be an investigation to understand why.

2cfe09b

do you want me to document that we have done a quick investigation and ran into a dead-end? I don't think it's worth investigating at this time.

[TonyXiang8787]

We don't need to mention anything in the docs. This is minor implementation choice.

We can later in the future to further benchmark why EIGEN is slower.

[TonyXiang8787]

This is for me an unusual way of describing the problem? Usually ,we are talking about 9 blocks in the process. But you only mentioned 4.

[
[a11, a12, a13],
[a21, a22, a23],
[a31, a32, a33],
]

Where
a11: block containing handled L/U combination where the pivots have been already handled
a12: already handled U part just above the current pivot
a13: already handled U part above the rest pivots to be handled
a21: already handled L part just in the left of current pivot
a22: current pivot
a23: current U part
a31: already handled L part in the left of the rest pivots to be handled
a32: current L part
a33: rest of the matrix which will be adjusted by current elimination process

[mgovers]

what do you mean i mentioned 4 blocks? i described the process for arbitrary block size. sym/asym is one thing, but the specific solvers can also add more dimensions to the block size. E.g. for NRPF, we treat the P and Q components separately, so it's not 3x3 but 9x9, right?

Since the process is iterative, you can step in at any pivot element $p$ and treat the rest as described until the next iteration $p$ needs to start.

[TonyXiang8787]

I mean, we should not omit the part which have been already factorized. We have to show the entire matrix in place, including the upper and left part of the matrix which already have the factorized L/U values.

[TonyXiang8787]

When you at pivot $k$, we should also show the rows and columns of $0 \to k-1$.

[TonyXiang8787]

1. Why do you use underscore? Underscore in EE means specifically complex numbers.
2. Same comment as above. Using 4 blocks instead of 9 blocks is confusing for the readers. Reader might think this is only applicable for the first pivot.

[mgovers]

I think i now understand what you mean. I'm describing the process at some point in the iteration during treatment of pivot element $p$.

I can indeed extend it, but the $p-1$ row and column do not contribute to this part of the process. What do you recommend?

[TonyXiang8787]

The $0 \to p-1$ does contribute in the process. If you permute columns/rows to get the max pivot for element $p$. The previous columns and rows have to be permuted as well.

Also, for the reader, include previous rows/columns shows that the factorization is in the middle of the process instead of just the beginning.

[TonyXiang8787]

This is not correct and was actually one of the attempts in the prototyping. It did not work and produced the same numerical instability as doing the explicit inverse.

The real implementation now, is we do not fully solve $x_{k;p}$, but only apply the two separate triangular solve in the current L and U part.

We apply lower-triangular solve in the U part which is right to the current pivot:

power-grid-model/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp

Lines 313 to 329 in 8d410ae

	// for block matrix
	// calculate U blocks in the right of the pivot, in-place
	// L_pivot * U_pivot,k = P_pivot * A_pivot,k k > pivot
	if constexpr (is_block) {
	for (Idx u_idx = pivot_idx + 1; u_idx < row_indptr[pivot_row_col + 1]; ++u_idx) {
	Tensor& u = lu_matrix[u_idx];
	// permutation
	u = (block_perm.p * u.matrix()).array();
	// forward substitution, per row in u
	for (Idx block_row = 0; block_row < block_size; ++block_row) {
	for (Idx block_col = 0; block_col < block_row; ++block_col) {
	// forward substract
	u.row(block_row) -= pivot(block_row, block_col) * u.row(block_col);
	}
	}
	}
	}

We then apply upper-triangular solver in the L part which is below the current pivot:

power-grid-model/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp

Lines 342 to 366 in 8d410ae

	if constexpr (is_block) {
	// for block matrix
	// calculate L blocks below the pivot, in-place
	// L_k,pivot * U_pivot = A_k_pivot * Q_pivot k > pivot
	Tensor& l = lu_matrix[l_idx];
	// permutation
	l = (l.matrix() * block_perm.q).array();
	// forward substitution, per column in l
	// l0 = [l00, l10]^T
	// l1 = [l01, l11]^T
	// l = [l0, l1]
	// a = [a0, a1]
	// u = [[u00, u01]
	// [0 , u11]]
	// l * u = a
	// l0 * u00 = a0
	// l0 * u01 + l1 * u11 = a1
	for (Idx block_col = 0; block_col < block_size; ++block_col) {
	for (Idx block_row = 0; block_row < block_col; ++block_row) {
	l.col(block_col) -= pivot(block_row, block_col) * l.col(block_row);
	}
	// divide diagonal
	l.col(block_col) = l.col(block_col) / pivot(block_col, block_col);
	}
	} else {

This is happening in the whole LU factorization part.

In the solve part, we do the rest lower/upper triangular solve on the actual b and y.

[sonarqubecloud]

Quality Gate passed

Issues
0 New issues
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
0.0% Duplication on New Code

See analysis details on SonarQube Cloud