Add Q-less tutorial

[MaxenceGollier]

Added two examples for solving
$AA^T y = b$ and $A^T A x = A^Tb$ without $Q$.

What do you think @amontoison ?

[dpo]

Please add text to describe that the examples are showing. It would also be useful to show how to use iterative refinement and compare the residuals and errors with and without it.

[amontoison]

Thanks @MaxenceGollier for adding a tutorial!

I have the following suggestions:

• Use m and n instead of 7 and 5;
• Use zeros instead of similar so that the user knows the shape of the vectors that should be allocated to solve the problems. We might want to do this before that we compute the right-hand side;
• Specify the option of not storing Q before the analysis so that QRMumps doesn't allocate memory to store the Q factor. You probably just don't update the Q factor if you do it after the analysis and before the factorization, but by then it's too late.

For the iterative refinement, we need to use get_factors to dump the sparse factor in Julia and perform the sparse triangular solves directly in Julia.
We should do that in another tutorial, but it will be easier if Alfredo (@abuttari) supports it directly in qr_mumps and its C interface.

[abuttari]

Thanks @MaxenceGollier !
I believe there is a mistake in the first example: y_star and y\_2 have to be of size n(=7) instead of m(=5). In fact you cannot compute norm(A*y\_2-b) because of incompatible dimensions. The actual least-norm solution is y\_3=A'*y\_2.
Have I missed something?

[abuttari]

For the iterative refinement, we need to use get_factors to dump the sparse factor in Julia and perform the sparse triangular solves directly in Julia. We should do that in another tutorial, but it will be easier if Alfredo (@abuttari) supports it directly in qr_mumps and its C interface.

@amontoison why can't you do the solves using qrm_solve?

[MaxenceGollier]

Thanks @MaxenceGollier ! I believe there is a mistake in the first example: y_star and y\_2 have to be of size n(=7) instead of m(=5). In fact you cannot compute norm(A*y\_2-b) because of incompatible dimensions. The actual least-norm solution is y\_3=A'*y\_2. Have I missed something?

In the first example, $A \in \mathbb{R}^{5 \times 7}$ and then, $AA^T \in \mathbb{R}^{5 \times 5}$. Hence, $b, y^*, y_2$ all lie in $\mathbb{R}^5$

[abuttari]

Thanks @MaxenceGollier ! I believe there is a mistake in the first example: y_star and y\_2 have to be of size n(=7) instead of m(=5). In fact you cannot compute norm(A*y\_2-b) because of incompatible dimensions. The actual least-norm solution is y\_3=A'*y\_2. Have I missed something?

In the first example, A∈R5×7 and then, AAT∈R5×5. Hence, b,y∗,y2 all lie in R5

@MaxenceGollier I get that. What I'm saying is that $y_2$ cannot be the solution of the problem because the solution has to be in $\mathbb{R}^7$

[amontoison]

For the iterative refinement, we need to use get_factors to dump the sparse factor in Julia and perform the sparse triangular solves directly in Julia. We should do that in another tutorial, but it will be easier if Alfredo (@abuttari) supports it directly in qr_mumps and its C interface.

@amontoison why can't you do the solves using qrm_solve?

You're right, we can. I was thinking of something else when I answered. We also don't need to call get_factors.

[MaxenceGollier]

What do you think now @dpo @amontoison @abuttari ?
$$y^*$$ is now in $$\mathbb{R}^7$$ as suggested.

[dpo]

Let's illustrate iterative refinement here, or else this tutorial is misleading. Also add a reference to Björck's paper where it is justified.

[dpo]

@MaxenceGollier

[MaxenceGollier]

@dpo I have added iterative refinement now.

[dpo]

I added suggestions to make the examples look more like a tutorial and be more didactic. Also, there seem to be extra steps.

Please modify example 2 accordingly.

[dpo]

I’m confused. This step doesn’t belong in the solution of the least-norm problem.

[MaxenceGollier]

See this comment. I am not sure if this answers your question @dpo ?

Thanks @MaxenceGollier ! I believe there is a mistake in the first example: y_star and y\_2 have to be of size n(=7) instead of m(=5). In fact you cannot compute norm(A*y\_2-b) because of incompatible dimensions. The actual least-norm solution is y\_3=A'*y\_2. Have I missed something?

[dpo]

Sorry, this needs to be adjusted. The method of seminormal equations for the least-norm problem is analyzed in by Paige in https://www.ams.org/journals/mcom/1973-27-122/S0025-5718-1973-0331745-1.

For normal equations, the stability with and without Q is similar.

But it's not the same for least squares. It's only for least squares that we need iterative refinement. I think it would be useful to write functions that solve those problems and hide the details.

[MaxenceGollier]

I will add iterative refinement in a separate PR.

[dpo]

I'm saying the text needs to be adjusted too.

[MaxenceGollier]

I think this is ready for a merge now, what do you think @dpo ?