# Solve Large Scale Underdetermined Linear Equation with per Element Equality Constraint

**URL:** <https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655>\
**Category:** Numerics\
**Tags:** linearalgebra, optimization, sparse, convex-optimization, linearsolve\
**Created:** [September 20, 2024, 7:01pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655 "2024-09-20T19:01:03Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 20, 2024, 7:01pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/1 "2024-09-20T19:01:03Z")

</div>

I have the following linear system to solve:

\boldsymbol{A} \boldsymbol{x} = \boldsymbol{0}, \quad \text{subject to} \; {x}\_{i} = {v}\_{i} \; \forall i \in \mathcal{V}

Where \boldsymbol{A} \in \mathbb{R}^{m \times n} with m \< n.  
The matrix \boldsymbol{A} is a banded matrix (5-9 bands).  
The set \mathcal{V} is a sub set of indices to force equality on a sub set of values of \boldsymbol{x}.

I wonder how to solve this in the most efficient way.

My current approach is to make this a Least Squares Problem:

\begin{align} \arg \min\_{\boldsymbol{x}} \quad & \frac{1}{2} \boldsymbol{x}^{T} \boldsymbol{B} \boldsymbol{x} \\ \text{subject to} \quad & \begin{aligned} {x}\_{i} & = {v}\_{i}, \; \forall i \in \mathcal{V} \\ \end{aligned} \end{align}

Where \boldsymbol{B} = \boldsymbol{A}^{T} \boldsymbol{A} is a Symmetric Positive Semi Definite (SPSD) matrix.

By defining a permutation matrix \boldsymbol{P} one can arrange the elements:

\begin{bmatrix} \boldsymbol{x}\_{\mathcal{U}} \\ \boldsymbol{x}\_{\mathcal{V}} \end{bmatrix} = \boldsymbol{P} \boldsymbol{x}

Which yields a partitioned formulation of the problem:

\boldsymbol{x}^{T} \boldsymbol{B} \boldsymbol{x} = \boldsymbol{x}^{T} \boldsymbol{P}^{T} \left( \boldsymbol{P} \boldsymbol{B} \boldsymbol{P}^{T} \right) \boldsymbol{P} \boldsymbol{x} = \begin{bmatrix} \boldsymbol{x}\_{\mathcal{U}}^{T} & \boldsymbol{x}\_{\mathcal{V}}^{T} \end{bmatrix} \begin{bmatrix} \boldsymbol{B}\_{\mathcal{U}} & \boldsymbol{C} \\ \boldsymbol{C} & \boldsymbol{B}\_{\mathcal{V}} \end{bmatrix} \begin{bmatrix} \boldsymbol{x}\_{\mathcal{U}} \\ \boldsymbol{x}\_{\mathcal{V}} \end{bmatrix}

Since only \boldsymbol{x}\_{\mathcal{U}} is sought, the problem becomes:

\arg \min\_{\boldsymbol{x}\_{\mathcal{U}}} \frac{1}{2} \boldsymbol{x}\_{\mathcal{U}}^{T} \boldsymbol{B}\_{\mathcal{U}} \boldsymbol{x}\_{\mathcal{U}} + \boldsymbol{x}\_{\mathcal{v}}^{T} \boldsymbol{C} \boldsymbol{x}\_{\mathcal{U}}

Which is a convex, smooth and constraint free problem.  
Its minimization is given where the gradient vanishes:

\boldsymbol{B}\_{\mathcal{U}} \hat{\boldsymbol{x}}\_{\mathcal{U}} = - \boldsymbol{C} \boldsymbol{x}\_{\mathcal{V}}

Which a simple linear system with SPSD matrix.

The advantages of this solution:

1. Smaller constrained free system.
2. SPSD system.

The weakness is being built by \boldsymbol{A}^{T} \boldsymbol{A} which increase the condition number.

Is there a better approach?  
Usually I’d solve the last equation using the Conjugate Gradient solver in [`Krylov.jl`](https://github.com/JuliaSmoothOptimizers/Krylov.jl). Is there an optimized solver for this?

**Remark** : I also [posted the question on StackExchange Computational Science](https://scicomp.stackexchange.com/questions/44550).

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 20, 2024, 7:25pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/2 "2024-09-20T19:25:58Z")

</div>

Why can’t you just delete `A[:, V]` and `x[V]` and solve the smaller system?

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [September 20, 2024, 8:45pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/3 "2024-09-20T20:45:48Z")

</div>

@amontoison this one is for you

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [September 20, 2024, 9:26pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/4 "2024-09-20T21:26:11Z")

</div>

Hi @RoyiAvital !  
Is your problem not equivalent to solve something similar to  
[A; I] [x; y] = [0; v] ?

The block I should be simplified.

Let call E the matrix needed for your problem where some rows of I are removed.  
If [A; E] has more rows than columns then use `lsmr` ohherwise use `craigmr`.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 21, 2024, 5:37am UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/5 "2024-09-21T05:37:36Z")

</div>

> [@Oscar\_Smith](#):
>
> Why can’t you just delete `A[:, V]` and `x[V]` and solve the smaller system?

Their value is used for the value of other elements.  
If I got your idea right, their value won’t cripple into other elements.

> [@amontoison](#):
>
> Hi @RoyiAvital !  
> Is your problem not equivalent to solve something similar to  
> [A; I] [x; y] = [0; v] ?
> 
> The block I should be simplified.
> 
> Let call E the matrix needed for your problem where some rows of I are removed.  
> If [A; E] has more rows than columns then use `lsmr` ohherwise use `craigmr`

Indeed, the original problem is given by:

\begin{bmatrix} \boldsymbol{A} \\ \boldsymbol{E} \end{bmatrix} \begin{bmatrix} \boldsymbol{x}\_{\mathcal{U}} \\ \boldsymbol{x}\_{\mathcal{V}} \end{bmatrix} = \begin{bmatrix} 0 \\ \boldsymbol{x}\_{\mathcal{V}} \end{bmatrix}

Where:

- \boldsymbol{E} = \begin{bmatrix} \boldsymbol{0} & \boldsymbol{I} \end{bmatrix}.
- \mathcal{U} = \left\{ 1, 2, \ldots, n \right\} \setminus \mathcal{V}.

But it is a larger system (n \times n) which is neither symmetric nor PSD. I thought it would be slower to solve.  
Does the solvers will be able to take advantage of its structure?  
Maybe I need to build a specialized operator?

I just thought in literature there are optimized solvers or approaches for such structure.

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [September 21, 2024, 6:05am UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/6 "2024-09-21T06:05:39Z")

</div>

@RoyiAvital It will not be slower because optimized Krylov solvers have been developped for rectangular problems (least-squares and least-norm).

`lsmr` and `craigmr` are equivalent to `minres` (in exact arithmetic) on the normal equations `AA' x = A'b` and `AA'y = b`. The optimality conditions of the least-squares and least-norm problems.

For information, `lsqr` and `craig` and equivalent to `cg` on the normal equations above.

Note that it’s way more stable than using `cg` and `minres` directly on the normal equations.  
It’s not recommended to use symmetric solvers on normal equations directly.

You create an operator if you don’t want to store [A; E] and just the blocks A and E. However, you will not gain performance, just storage.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 21, 2024, 7:25am UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/7 "2024-09-21T07:25:54Z")

</div>

OK, I will give it a try.  
One delicate thing (At least when using _direct solvers_), The LS approach will always result in a solution with the equality forced. The solution with \boldsymbol{E}, in case no solution in the range of the model matrix will result in LS solution. So the equality will be approximated.

I will create a short script and share results.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 21, 2024, 1:39pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/8 "2024-09-21T13:39:02Z")

</div>

I implemented both methods from above.  
The matrix is based on the definition in [Build the Affinity Matrix of Edge Preserving Multiscale Image Decomposition based on Local Extrema](https://dsp.stackexchange.com/questions/95071).  
The image size was `250 x 250`. The size of of the set \left| \mathcal{V} \right| = 13747.  
The results I got, using Julia:

- Direct Solvers:
  - Direct Solution (Solving the formulation with \boldsymbol{A}):
    - Error: 0 (Will be used as reference).
    - Run Time: `0.14130` [Sec].

  - LS Solution (Solving the formulation with \boldsymbol{B})
    - Error: 4e-12.
    - Run Time: `0.21624` [Sec].

- Iterative Solvers:
  - LS Solution based on Conjugate Gradient:
    - Error (Max Absolute vs. Direct Solution): `0.000047055766`.
    - Run Time: `8.86836` [Sec].

  - Direct Solution based on LSMR:
    - Error (Max Absolute vs. Direct Solution): `0.086388343173`.
    - Run Time: `4.27441` [Sec].

  - Direct Solution based on LSQR:
    - Error (Max Absolute vs. Direct Solution): `0.001536492862`.
    - Run Time: `7.62469` [Sec].

  - Direct Solution based on CRAIG:
    - Error (Max Absolute vs. Direct Solution): `0.000066854811`.
    - Run Time: `9.38470` [Sec].

  - Direct Solution based on CGLS:
    - Error (Max Absolute vs. Direct Solution): `0.000083813606`.
    - Run Time: `8.60961` [Sec].

From my point of view, if direct solver is feasible, the direct formulation will be the fastest.  
If iterative method is needed, probably go with `CG` if accuracy is required or `LSMR` if speed is required.

It might be different with more tweaking, but it seems they all are in the same ballpark.

All iterative solvers are form [`Krylov.jl`](https://github.com/JuliaSmoothOptimizers/Krylov.jl) with the default settings.  
Run time measured by [`BenchmarkTools.jl`](https://github.com/JuliaCI/BenchmarkTools.jl) using `@belapsed`.

**Remark** : It seems most of the time in the LS method is the calculation of \boldsymbol{B} = \boldsymbol{A}^{T} \boldsymbol{A}. Once there is a threaded function for sparse matrix - sparse matrix product it will be faster than the direct method. A PoC using [`SuiteSparseGraphBLAS.jl`](https://github.com/JuliaSparse/SuiteSparseGraphBLAS.jl) showed it is faster. Yet in my case the matrices are [`SparseMatrixCSC{Float64, Int32}` which is not supported in `SuiteSparseGraphBLAS.jl`](https://github.com/JuliaSparse/SuiteSparseGraphBLAS.jl/issues/144). I [asked for it in `MKLSparse.jl` a s well](https://github.com/JuliaSparse/MKLSparse.jl/issues/53).

---

<div class="post-metadata">

**Author:** ![amontoison](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/amontoison/32/218741_2.png) [@amontoison](https://discourse.julialang.org/u/amontoison)\
**Post date:** [September 21, 2024, 3:57pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/9 "2024-09-21T15:57:17Z")

</div>

A’A is only computed if you use `cg` on the normal equations, otherwise you should provide `A` as input.  
Did you try with `using MKL` and `using MKLSparse` before `using Krylov`?

I recently updated `MKLSparse` such that sparse matrix × dense vector products are multithreaded for the both types of integer (Int32 and Int64).

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 21, 2024, 4:53pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/10 "2024-09-21T16:53:40Z")

</div>

I did provide only \boldsymbol{A}.  
I tested all of them using the same environment, so it is fair comparison as `MKLSparse.jl` will benefit all of them (Iterative methods).

I don’t think it will change the conclusions:

1. Using Direct Solver: If there is a fast Sparse Matrix - Sparse Matrix multiplication, use the LS (Normal Equations). Otherwise the direct form of the problem.
2. Using Iterative Solver: For accuracy, use `cg()`. For speed use `lsmr()`.

Does it surprise you? Want me to check something more?

The code is available on my [StackExchange Computational Science GitHub Repository](https://github.com/RoyiAvital/StackExchangeCodes) (Look at the `ComputationalScience\Q44552` folder).

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [September 23, 2024, 11:30am UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/11 "2024-09-23T11:30:27Z")

</div>

I might be taking this a step back after you have already received some more specific responses, but I’m not sure I fully understood the problem. Adding some further partitioning, you essentially state that your large system is of the form

\begin{bmatrix} \boldsymbol{A}\_{\mathcal{U}} & \boldsymbol{A}\_{\mathcal{V}} \\ 0 & \boldsymbol{I} \end{bmatrix} \begin{bmatrix} \boldsymbol{x}\_{\mathcal{U}} \\ \boldsymbol{x}\_{\mathcal{V}} \end{bmatrix} = \begin{bmatrix} 0 \\ \boldsymbol{x}\_{\mathcal{V}} \end{bmatrix},

where the matrix is n\times n. This implies that I is (n-m)\times (n-m) and that \boldsymbol{A}\_{\mathcal{U}} is m\times m. If \boldsymbol{A}\_{\mathcal{U}} has rank m, then you can set \boldsymbol{x}\_{\mathcal{U}} = \boldsymbol{v} for some prescribed \boldsymbol{v} as in your original post. Then you have

\boldsymbol{A}\_{\mathcal{U}} \boldsymbol{x}\_{\mathcal{U}} = -\boldsymbol{A}\_{\mathcal{V}} \boldsymbol{v}

which would be a square nonsingular system. Since \boldsymbol{B}\_{\mathcal{U}} = \boldsymbol{A}\_{\mathcal{U}}^T \boldsymbol{A}\_{\mathcal{U}}, your other formulation would also be positive definite, although I don’t see much reason to form \boldsymbol{B}\_{\mathcal{U}} when you can solve the system with \boldsymbol{A}\_{\mathcal{U}} directly.

On the other hand, if \boldsymbol{A}\_{\mathcal{U}} has rank less than m, either A has rank less than m or you made an unfortunate choice of columns to include in \boldsymbol{A}\_{\mathcal{U}}. (I’m not sure if you have any choice here…) So assuming \boldsymbol{A}\_{\mathcal{U}} is singular, I am wondering which case it is: Does A have rank less than m or are you just stuck with a bad choice of columns in \boldsymbol{A}\_{\mathcal{U}} that you are not free to change? For singular \boldsymbol{A}\_{\mathcal{U}} is there some reason to believe a solution exists for your particular choice of \boldsymbol{v}? If it exists are you looking for a minimum norm solution or any solution? If it doesn’t are you looking for a minimum norm least squares solution?

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [September 23, 2024, 6:11pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/12 "2024-09-23T18:11:00Z")

</div>

> [@Oscar\_Smith](#):
>
> Why can’t you just delete `A[:, V]` and `x[V]` and solve the smaller system?

@RoyiAvital this seems to me like the most sensible suggestion here. Delete columns of A corresponding to variables in V and solve a thinner least-squares problem. Why would you iterate on fixed variables? Either use LSMR or Golub/Riley on the eliminated problem.

ps: it doesn’t help to spread your questions across several threads. I did not see this before.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 24, 2024, 6:07am UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/13 "2024-09-24T06:07:17Z")

</div>

> [@mstewart](#):
>
> I might be taking this a step back after you have already received some more specific responses, but I’m not sure I fully understood the problem. Adding some further partitioning, you essentially state that your large system is of the form

This is actually how I solve it in [my code for Edge Preserving Multiscale Image Decomposition based on Local Extrema](https://github.com/RoyiAvital/StackExchangeCodes/blob/f46513fa2836615032be000da07686629f8a51e1/SignalProcessing/Q95071/Q95071.jl#L151).  
In practice in some case those Image Processing problems (Form the _pre historic era_) you get SPSD matrix to begin with (Like \boldsymbol{B}, The Laplacian of the graph). Yet in this formulation the matrix is not even symmetric. Hence “Pseudo Laplacian”.

I like implementing those papers in my free time.  
They are probably useless these days, yet they make the thinking wheels move.

> [@mstewart](#):
>
> On the other hand, if \boldsymbol{A}_{\mathcal{U}}AU\boldsymbol{A}_{\mathcal{U}} has rank less than mmm, either AAA has rank less than mmm or you made an unfortunate choice of columns to include in \boldsymbol{A}_{\mathcal{U}}AU\boldsymbol{A}_{\mathcal{U}}. (I’m not sure if you have any choice here…) So assuming \boldsymbol{A}_{\mathcal{U}}AU\boldsymbol{A}_{\mathcal{U}} is singular, I am wondering which case it is: Does AAA have rank less than mmm or are you just stuck with a bad choice of columns in \boldsymbol{A}_{\mathcal{U}}AU\boldsymbol{A}_{\mathcal{U}} that you are not free to change? For singular \boldsymbol{A}_{\mathcal{U}}AU\boldsymbol{A}_{\mathcal{U}} is there some reason to believe a solution exists for your particular choice of \boldsymbol{v}v\boldsymbol{v}? If it exists are you looking for a minimum norm solution or any solution? If it doesn’t are you looking for a minimum norm least squares solution?

The formulation of the matrix means it is rank deficient.  
As rows of \mathcal{U} must have zero sum.  
Hence any vector which has constant values is an eigen vector of \boldsymbol{A}\_{\mathcal{U}} with eigen value of 0.

> [@dpo](#):
>
> @RoyiAvital this seems to me like the most sensible suggestion here. Delete columns of A corresponding to variables in V and solve a thinner least-squares problem

I might overlooked this. But unless you do the formulation done for \boldsymbol{B} (Spelled out in @mstewart 's post) I don’t see how deleting can work. You must build a different set of equations. Deleting won’t work as the values are needed.

> [@dpo](#):
>
> ps: it doesn’t help to spread your questions across several threads. I did not see this before.

The other thread is solely on applying _Cholesky Decomposition_ to a rank deficient sparse SPSD matrix. I was wondering how to handle it given pivoting is not available for the sparse case.  
It drifted to a similar discussion as people, rightfully, asked for the context.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [September 24, 2024, 8:53am UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/14 "2024-09-24T08:53:53Z")

</div>

> [@RoyiAvital](#):
>
> The formulation of the matrix means it is rank deficient. As rows of \mathcal{U} must have zero sum. Hence any vector which has constant values is an eigen vector of \boldsymbol{A}\_{\mathcal{U}} with eigenvalue of 0.

That is a good explanation. It means your solution won’t be unique. But is the system expected to be consistent? What sort of solution are you looking for? A minimum norm least squares solution?

> [@RoyiAvital](#):
>
> I might overlooked this. But unless you do the formulation done for \boldsymbol{B} (Spelled out in @mstewart 's post) I don’t see how deleting can work. You must build a different set of equations. Deleting won’t work as the values are needed.

I think this was probably the same suggestion as mine and I was just a bit more explicit about the details. After I made my post, I looked back and thought I should have pointed out that I wasn’t the first person making the suggestion. If I understand correctly now, you pick \boldsymbol{x}\_{\mathcal{V}} and solve for \boldsymbol{x}\_{\mathcal{U}} either way. So I’m still not seeing what benefit you get working with \boldsymbol{B} instead of working with \boldsymbol{A}\_{\mathcal{U}}.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 24, 2024, 9:29am UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/15 "2024-09-24T09:29:58Z")

</div>

> [@mstewart](#):
>
> I think this was probably the same suggestion as mine and I was just a bit more explicit about the details. After I made my post, I looked back and thought I should have pointed out that I wasn’t the first person making the suggestion. If I understand correctly now, you pick \boldsymbol{x}_{\mathcal{V}}xV\boldsymbol{x}_{\mathcal{V}} and solve for \boldsymbol{x}_{\mathcal{U}}xU\boldsymbol{x}_{\mathcal{U}} either way. So I’m still not seeing what benefit you get working with \boldsymbol{B}B\boldsymbol{B} instead of working with \boldsymbol{A}_{\mathcal{U}}AU\boldsymbol{A}_{\mathcal{U}}.

When I did it first I just went with intuition.  
Later I just proved it mathematically using permutation matrices (Sub set of the rows of a permutation matrix).  
I overlooked on just employing the simple choice and creating the equation as you did.  
Indeed, I now think this is what @Oscar_Smith meant. Yet I took it literally as removing the rows.

I guess the only benefit comes to the simple balance when dealing with the LS System.  
You may gain some speed in the case m \< n solving the normal equations but you pay with accuracy as you square the condition number.  
In some cases it worth doing (Classic application is Harris Corner Detection and Optical Flow), in others it might be too dangerous. I guess in real world it depends on the given time budget.

---

<div class="post-metadata">

**Author:** ![dpo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpo/32/3335_2.png) [@dpo](https://discourse.julialang.org/u/dpo)\
**Post date:** [September 24, 2024, 6:37pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/16 "2024-09-24T18:37:09Z")

</div>

> [@RoyiAvital](#):
>
> I might overlooked this. But unless you do the formulation done for \boldsymbol{B}B\boldsymbol{B} (Spelled out in @mstewart 's post) I don’t see how deleting can work. You must build a different set of equations. Deleting won’t work as the values are needed.

If V is the set of fixed variables (that’s what they are called in optimization) and F is the set of remaining (free) variables, partition WLOG,

A = \begin{bmatrix} A\_F & A\_V \end{bmatrix}

and x accordingly. Then,

Ax = 0 \Longleftrightarrow A\_F x\_F = -A\_V x\_V.

Thus, if A\_F is still underdetermined, you can minimize \| x\_F \| subject to the eliminated system above.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [September 25, 2024, 5:10pm UTC](https://discourse.julialang.org/t/solve-large-scale-underdetermined-linear-equation-with-per-element-equality-constraint/119655/17 "2024-09-25T17:10:56Z")

</div>

Indeed.  
This is exactly the same trick done by me in the code and spelled by @mstewart .  
I didn’t understand @Oscar_Smith 's post that way, hence I overlooked it.