Linear Algebra

Diameter theorem for polytopes and polyhedra: Every diameter of an \(n\)-dimensional polytope or polyhedron given by \(m\) constraints with \(m,n\in{}^{\omega}\mathbb{N}_{\ge 2}\) is at most \(2(m + n – 3)\).

Proof: At most \(\acute{m}\) hyperplanes can be combined to form an incomplete cycle of dimension \(2\), and there are at most \(n-2\) options to deviate sideways in the remaining dimensions. Traversing each minimal segment requires at most two edges and thus yields the factor \(2\).\(\square\)

Theorem on Strassen’s algorithm: For sufficiently large \(n := 2^{\ell}\), \(\ell \in {}^{\nu}\mathbb{N}^*\), \(ß := {}_27\) and \(A \in {}^{\nu}\mathbb{C}^{n \times n}\), the GS proves the reduction of the running time \(T(n) = \mathcal{O}(n^{ß})\)¹ by approximately \(\tilde{3}\) for the computation of \[ AA^H = \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix} \begin{pmatrix} A_{11}^H & A_{21}^H \\ A_{12}^H & A_{22}^H \end{pmatrix} = \begin{pmatrix} A_{11}A_{11}^H+A_{12}A_{12}^H & A_{11}A_{21}^H+A_{12}A_{22}^H \\ A_{21}A_{11}^H+A_{22}A_{12}^H & A_{21}A_{21}^H+A_{22}A_{22}^H \end{pmatrix}.\;\square \] Theorem on fast matrix multiplication: With \(A\) and \(B = (b_{ij}) \in {}^{\nu}\mathbb{C}^{n \times n}\) as above and \(n\) as before, \(AB\) can be computed for \[ s := \min\{2^k > \max_{i,j} (|a_{ij}|^2, |b_{ij}|^2)\hat{n} : k \in {}^{\nu}\mathbb{Z}\} \] from \[ A_{11}(s B_{11} + B_{12}),\; A_{21}(s B_{11} + B_{12}),\; A_{12}(s B_{21} + B_{22}),\; A_{22}(s B_{21} + B_{22}) \] by computing modulo \(s\) in running time \(\mathcal{O}(p^2)\) with \(p = \ell n\) (see Matrix Multiplication).\(\square\)

Revised INLIF and Farkas certificates

The full INLIF tableau makes all current coefficients explicitly visible. This is advantageous for the derivation, for traceability and for the tableau trace. For large or sparse systems, however, it is preferable not to store the entire transformed tableau, but only the current basis structure and the rows and columns needed for the next step. This representation is called \emph{revised INLIF} in what follows.

The starting point is the system \(Ax\le b\). Introducing slack variables \(s\ge0\) gives \(Ax+s=b\). With \(H=(A\; I)\) and \(q=(x,s)^\top\), the extended system of equations is \(Hq=b\). In free INLIF, the original variables \(x_j\) are initially free variables, whereas the slack variables \(s_i\) are cone variables. At each point in time, \(m\) variables stand on the row side. The corresponding columns of \(H\) form the current basis matrix \(\mathcal B\). The remaining columns form the current column matrix \(\mathcal N\). Thus \( \mathcal B q_R+\mathcal N q_C=b \) and consequently \( q_R=\mathcal B^{-1}b-\mathcal B^{-1}\mathcal N q_C. \) Setting \(\beta:=\mathcal B^{-1}b\) and \(G:=\mathcal B^{-1}\mathcal N\), the full tableau is \( q_R=\beta-Gq_C. \)

Revised INLIF does not form the matrix \(G\) in full. Instead, only \(\mathcal B\), or a factorisation or product-form representation of \(\mathcal B\), is stored. The current right-hand side is obtained from \(\mathcal B\beta=b\). An active test row is violated if the corresponding row variable is a cone variable and \(\beta_k<0\). Row variables that have arisen from originally free variables are, by contrast, pure back-substitution rows and are not tested for non-negativity.

For pivot selection in a violated row \(k\), the whole tableau row is not computed. Instead, the transposed basis system \(\mathcal B^\top \lambda_k=e_k\) is solved. Then \( \lambda_k^\top=e_k^\top \mathcal B^{-1}. \) For a current column variable with original column \(h_j\) in \(H\), the corresponding tableau coefficient is obtained from \( G_{kj} = e_k^\top\mathcal B^{-1}h_j = \lambda_k^\top h_j. \) This is the central revised-INLIF formula: \[ \boxed{G_{kj}=\lambda_k^\top h_j.} \] A free column variable is admissible as a pivot if \(G_{kj}\ne0\). A cone column variable is admissible as a pivot only if \(G_{kj}<0\). Thus the pivot logic of the full INLIF tableau remains unchanged; only the computation of the coefficients is implicit.

Once a pivot column \(p\) has been chosen, the transformed column is computed from \(\mathcal B u=h_p\). Then \(u=\mathcal B^{-1}h_p\) and, in particular, \(u_k=G_{kp}\). The new basis is obtained by having the column variable \(q_{C_p}\) replace the previous row variable \(q_{R_k}\). Writing \(\mathcal B’=\mathcal B E\), where \(E\) is the identity matrix whose \(k\)-th column has been replaced by \(u\), it follows that \( \mathcal B’^{-1}=E^{-1}\mathcal B^{-1}. \) The right-hand side can therefore be updated without a full tableau: \( \beta’=E^{-1}\beta. \) For \(i\ne k\), the structure of \(E\) explicitly yields \(\beta’_k=\beta_k/u_k\) and \( \beta’_i = \beta_i-u_i\beta’_k. \) Thus an eta update replaces the full tableau update.

The advantage does not lie in building a second dual system. Revised INLIF does not require \(A\) and \(A^\top\) to be maintained separately as two independent inequality systems. The transposed equations \(\mathcal B^\top\lambda_k=e_k\) and, in the LP case, \(\mathcal B^\top y=\breve c_R\), refer to the same current basis \(\mathcal B\). The dual information therefore arises from the cross-tableau structure itself and not from a separate INLIF treatment of the dual problem.

For an LP \[ \max\; e+c^\top x \qquad\text{subject to}\qquad Ax\le b,\quad x\ge0 \] the extended cost vector \(\breve c=(c,0)^\top\) is used. If \(\breve c_R\) are the costs of the current row variables, the objective multiplier \(y\) is determined from \( \mathcal B^\top y=\breve c_R. \) For a column variable \(q_j\) with column \(h_j\), the reduced coefficient is \( \overline c_j=\breve c_j-h_j^\top y. \) The dual inequalities are therefore not stored as a separate system, but read from the reduced coefficients of the current basis. For the original \(x\)-columns this gives \( \overline c_{x_j}=c_j-A_{\cdot j}^\top y, \) and for the original slack columns \(\overline c_{s_i}=-y_i\). The conditions \(\overline c_{x_j}\le0\) and \(\overline c_{s_i}\le0\) therefore correspond to \(A^\top y\ge c\) and \(y\ge0\). Together with a primal feasible \(x\) and \(b^\top y-c^\top x=0\), this gives a primal-dual optimality certificate.

The extraction of a Farkas certificate is especially natural in revised INLIF. Let \(k\) be a violated active test row with \(\beta_k<0\). For this row, the system \(\mathcal B^\top\lambda_k=e_k\) has already been solved for the pivot test. Hence \[ b^\top\lambda_k = \lambda_k^\top b = e_k^\top\mathcal B^{-1}b = \beta_k < 0. \] If the row is blocked, meaning that neither a free pivot nor an admissible cone pivot exists, then \(\lambda_k\) can be checked directly against the original system.

For originally free variables, a blocked row must satisfy \(A^\top\lambda_k=0\). Indeed, any non-zero coefficient of a free column variable would permit a free pivot. For cone variables, only non-negative coefficients may remain in the blocked row. Since the slack columns form the identity matrix, these coefficients are precisely the components of \(\lambda_k\). Thus \(\lambda_k\ge0\). Hence, if \[ \lambda_k\ge0, \qquad A^\top\lambda_k=0, \qquad b^\top\lambda_k<0 \] is verified, then \[ \boxed{ \lambda_k \text{ is a Farkas certificate for the infeasibility of } Ax\le b. } \]

Indeed, from a hypothetical \(x\) with \(Ax\le b\) one would obtain the contradiction \[ 0 = x^\top A^\top\lambda_k = \lambda_k^\top Ax \le \lambda_k^\top b < 0. \] For systems with the additional non-negativity condition \(x\ge0\), the corresponding Farkas check is \[ \lambda_k\ge0, \qquad A^\top\lambda_k\ge0, \qquad b^\top\lambda_k<0. \] For every hypothetical \(x\ge0\) with \(Ax\le b\), this would then give the contradiction \[ 0 \le x^\top A^\top\lambda_k = \lambda_k^\top Ax \le \lambda_k^\top b < 0. \]

Thus revised INLIF yields the following certified alternative: \[ \boxed{ \begin{array}{ll} \text{either} & \text{the current origin is feasible and yields an original point,}\\[1mm] \text{or} & \text{a blocked violated row yields a Farkas certificate.} \end{array} } \] In the LP case, a primal-dual optimality certificate is added as a third terminal case: \[ Ax\le b,\quad x\ge0,\quad A^\top y\ge c,\quad y\ge0,\quad b^\top y-c^\top x=0. \] Unboundedness, by contrast, is not inferred merely from an open improving column. It is accepted only if a ray \(d\) with \(d\ge0\), \(Ad\le0\) and \(c^\top d>0\) is verified in the original system.

The step bound follows from the monotonicity of the absorptions. Each free absorption removes a still free variable from the active residual system. Cone absorptions arise only from origin positions that have already been activated. Under the extremal rule and the non-reactivation of activated hyperplanes, the stable reduction classes therefore satisfy \[ N_{\mathrm{steps}} = N_{\mathrm{free}}+N_{\mathrm{cone}} \le n+n = \hat n. \] For LP s with \(n\) primal decision variables and \(m\) slack or row variables, the corresponding structural bound \(N_{\mathrm{absorption}}\le m+n\) is used, provided that every tableau variable is absorbed into the computed role at most once. Reaching this bound is not in itself a proof of infeasibility. Afterwards, either a primal-dual optimality certificate or a verified Farkas or unboundedness certificate must be present. Otherwise, the monotonicity or pivot invariant has been violated, and the case remains algorithmically open.

Let \(s\) denote the number of absorption steps and \(F_t\) the cost of a basis solve in step \(t\). If, in step \(t\), an average of \(q_t\) candidate columns is tested and the matrix has effective density \(\delta\), the working estimate for revised INLIF is \[ T_{\mathrm{rev}} = \mathcal O\!\left( L{\LARGE{\textbf{+}}}_{t=1}^{s} \left( F_t+q_t\delta m \right) \right). \] With full column scanning, \(q_t=\mathcal O(n)\). For \(m=\mathcal O(n)\) and \(s=\mathcal O(n)\), the scanning part is then \(\mathcal O(L\delta n^3)\). With local pivot testing and \(q_t=\mathcal O(1)\), this part formally drops to \(\mathcal O(L\delta n^2)\). The costs of the basis solves must additionally be taken into account. In the dense worst case, the familiar cubic order may remain. The decisive gain of revised INLIF therefore lies in reduced memory consumption, the use of sparsity and the fact that the dense tableau \(G=\mathcal B^{-1}\mathcal N\) does not have to be fully stored.\(\square\)

Problem statement and meaning of \(\operatorname{pro}\).
The problem statement is the same as for the INLIF method. The parameter \(\operatorname{pro}\in\{1,-1\}\) determines how the input data are interpreted.

For \(\operatorname{pro}=1\), the input problem is read as a maximising primal problem: \[ \max \; e+c^\top x \qquad\text{unter}\qquad Ax\le b,\qquad x\ge 0. \] For \(\operatorname{pro}=-1\), the internal primal-dual orientation is exchanged. The method then works with the transposed-negated tableau form, so that the internal objective function plays the same geometric role as in the maximising case. However, the mathematical objective value must not be output with the wrong sign. Therefore the internally scaled objective value is multiplied by \(\operatorname{pro}\) during back-transformation: \[ f_{\mathrm{orig}} = e+\operatorname{pro}\,\kappa f_{\mathrm{int}}, \qquad \kappa>0. \] Thus the factor \(\operatorname{pro}\) belongs to the back-transformation of the internal objective value, not to the primal and dual feasibility conditions themselves. The test quantities \[ Ax\le b,\qquad A^\top y\ge c,\qquad b^\top y-c^\top x\ge 0 \] are always evaluated in the respective back-transformed original orientation.

Theorem.
The hybrid Intex–Merit method determines a primal-dual solution \((x^o,y^o)\) for every solvable LP of the above form, provided that the arithmetic precision used, the relaxation sequence, and the final CO resolve the active structure of the problem. The Intex core works by Inter-/Extrapolations. For a matrix density \(d\in[0,1]\) and a suitable binary precision scale \({}_2(\tilde\alpha\rho)\), the geometric Intex phase has, in the unstructured case, the order \[ \mathcal O\!\left({}_2(\tilde\alpha\rho)^2\,dmn\right). \] If the linear algebra is carried out in structured form by FTD/FSR, the respective structured operator cost takes the place of \(dmn\). Let \(z:=m+n.\)

Intex relaxation.
Instead of enforcing directly \[ Ax\le b,\qquad A^\top y\ge c,\qquad b^\top y-c^\top x=0 \] the Intex method considers, for \(r\in[0,\rho]\), a family of relaxed primal-dual polytopes \[ P_r := \left\{ (x,y)^\top\in{}^{\nu}\mathbb R_{\ge 0}^{z} : b^\top y-c^\top x\le r,\; Ax-b\le r_{m\times 1},\; c-A^\top y\le r_{n\times 1}, \right\}. \] The initial radius is chosen so that the origin \(0\) is certainly contained in \(P_\rho\). One possible choice is \[ \rho := s\left|\min\{b_1,\dots,b_m,-c_1,\dots,-c_n\}\right|, \qquad s\in]1,2], \] possibly after prior scaling and shifting of the data.

The original LP is solved exactly when a pair \((x,y)\in P_0\) is found. By strong duality², this is equivalent to \[ Ax\le b,\qquad x\ge0,\qquad A^\top y\ge c,\qquad y\ge0,\qquad b^\top y-c^\top x=0. \] Thus the same variables simultaneously solve \[ \max\{c^\top x:x\in{}^{\nu}\mathbb R_{\ge0}^n,\ Ax\le b\} \] and \[ \min\{b^\top y:y\in{}^{\nu}\mathbb R_{\ge0}^m,\ A^\top y\ge c\}. \]

Normalisation and scaling.
Before the geometric iteration, the objective function, right-hand sides, and constraint matrix are normalised. The aim is not to change the problem, but to obtain a balanced representation of the three types of error \[ Ax-b,\qquad c-A^\top y,\qquad b^\top y-c^\top x. \] The scaling is intended to prevent one of these three quantities from dominating numerically. To ensure that the objective function coupling in the gradient descent of the subsequent merit polishing neither disappears nor the hinge errors are numerically masked, an initial scalar linear transformation is applied. Specifically, the vector \(c\) and the constant offset \(e\) are scaled by the factor \(\lambda = \|c\|^{-1}\) (for \(c \neq 0\)).

Under this transformation, the primal solution \(x\) remains invariant, while the dual solution \(y\) and the absolute objective function values scale strictly linearly with \(\lambda\). At the end of the overall procedure, an exact back-transformation to the original physical space is performed. After the computation, all variables and objective values are transformed back to the original orientation. In particular, when \(\operatorname{pro} = -1\), the internally computed objective value is output with the correct sign.

Geometric Intex step.
Within a fixed polytope \(P_r\), an interior centre of mass is approximated. For each coordinate, a feasible interval is determined. For \(v \in P_r\), one approximately forms \[ v_k^* := \min\check v_k+\max\check v_k, \qquad k=1,\dots,z, \] This yields a geometric centre \(v=(x,y)^\top\), which remains inside the current relaxation and serves as a stable starting point for the next radius reduction.

The most recent centre movements are then extrapolated. For this purpose, a direction \(\Delta v\) is formed, and a maximally feasible step along \[ v(w)=v+w\Delta v,\qquad w\ge0, \] is determined. The new radius is obtained from the smallest relaxation still required: \[ r_{\mathrm{new}} = \min\{r\ge0:v(w)\in P_r\}. \] If \(r_{\mathrm{new}}=0,\) is achieved, a primal-dual optimal pair has been reached. If \(r_{\mathrm{new}}>0\) remains stably separated from zero, this indicates infeasibility, insufficient precision, or an unresolved active structure.

Intex as a feasible starting point.
The essential advantage of the Intex phase is that it does not search blindly in the exterior space. It works with relaxed feasible regions and typically provides a point which is already well oriented in primal-dual terms: \[ Ax-b\lesssim r_{m\times 1},\qquad c-A^\top y\lesssim r_{n\times 1},\qquad b^\top y-c^\top x\lesssim r. \] This point is not necessarily optimal yet, but it is often close to an active face of the optimal polytope. This is precisely where the Merit polishing phase starts.

Primal-dual Merit polishing and algebraic crossover

Residuals.
For \(x\ge0\) and \(y\ge0\), let \[ r^p(x):=\max(Ax-b,0)\in{}^{\nu}\mathbb R^m, \] \[ r^d(y):=\max(c-A^\top y,0)\in{}^{\nu}\mathbb R^n \] denote the primal and dual violations. The signed duality gap is \[ g(x,y):=b^\top y-c^\top x. \] For a maximising primal problem, only \[ g_+(x,y):=\max(g(x,y),0) \] is relevant as a positive gap violation. At the exact optimum, \[ r^p(x)=0,\qquad r^d(y)=0,\qquad g(x,y)=0. \] The unscaled control quantities are \[ p_\infty=\|r^p(x)\|_\infty, \qquad d_\infty=\|r^d(y)\|_\infty, \qquad g_{\mathrm{abs}}=|b^\top y-c^\top x|. \] The primal and dual objective values are \[ p^*=c^\top x, \qquad d^*=b^\top y. \] The constant offset \(e\) is added only at output.

Merit function.
The Merit function combines primal violation, dual violation, and duality gap: \[ \widehat{\Phi}(x,y) = \|r^p(x)\|_2^2 + \|r^d(y)\|_2^2 + \gamma^2\,g_+(x,y)^2, \qquad \gamma>0. \] In scaled variables \[ x=S_x\breve x,\qquad y=S_y\breve y \] and with diagonal residual weights \[ D_p=\operatorname{diag}(d_p), \qquad D_d=\operatorname{diag}(d_d) \] the scaled form is \[ \widehat\Phi(\breve x,\breve y) = \|D_p\max(AS_x\breve x-b,0)\|_2^2 + \|D_d\max(c-A^\top S_y\breve y,0)\|_2^2 + d_g^2 \max(b^\top S_y\breve y-c^\top S_x\breve x,0)^2. \] Optionally, a small regularisation term \[ \varepsilon\left(\|\breve x\|_2^2+\|\breve y\|_2^2\right) \] may be added.

Gradient.
With \[ w_i^p=d_{p,i}^2\,r_i^p(\breve x), \qquad w_j^d=d_{d,j}^2\,r_j^d(\breve y), \qquad w_g=d_g^2\,g_+(\breve x,\breve y) \] the gradients are \[ \frac{{\downarrow}\Phi}{{\downarrow}\breve x_j} = s_{x,j} \left( {\LARGE{\textbf{+}}}_{i=1}^m A_{ij}w_i^p – c_j w_g \right), \qquad j=1,\dots,n, \] \[ \frac{{\downarrow}\Phi}{{\downarrow}\breve y_i} = s_{y,i} \left( b_i w_g – {\LARGE{\textbf{+}}}_{j=1}^n A_{ij}w_j^d \right), \qquad i=1,\dots,m. \] The term \(b_iw_g\) is essential: the gradient with respect to \(y\) must contain the contribution of \(b^\top y\) in the duality gap. Without this contribution, the dual side is steered incorrectly.

Feasible step control.
Starting from the Intex point, a descent step is computed. In the simple projected variant, \[ \breve x’ = \Pi_{\mathbb R_+^n} \left( \breve x-\alpha\nabla_{\breve x}\widehat\Phi \right), \qquad \breve y’ = \Pi_{\mathbb R_+^m} \left( \breve y-\alpha\nabla_{\breve y}\widehat\Phi \right). \] In the feasible hybrid variant, a hard orthant projection is avoided during the Merit phase. Instead, a step is rejected or shortened if it would leave the feasible relaxation or violate non-negativity. This preserves the character of the Intex relaxation.

A step is accepted only if an Armijo-type condition is satisfied: \[ \widehat\Phi(\breve x’,\breve y’) \le \widehat\Phi(\breve x,\breve y) – c_1\alpha \left( \|\breve x’-\breve x\|_2^2+ \|\breve y’-\breve y\|_2^2 \right). \] Otherwise, \(\alpha\) is reduced.

Barzilai–Borwein step size.
After an accepted step, \[ s_k= \begin{pmatrix} \breve x^{\grave{k}}-\breve x^k\\[0.3em] \breve y^{\grave{k}}-\breve y^k \end{pmatrix}, \qquad q_k= \begin{pmatrix} \nabla_{\breve x}\widehat{\Phi}(\breve x^{\grave{k}},\breve y^{\grave{k}})-\nabla_{\breve x}\widehat{\Phi}(\breve x^k,\breve y^k)\\[0.3em] \nabla_{\breve y}\widehat{\Phi}(\breve x^{\grave{k}},\breve y^{\grave{k}})-\nabla_{\breve y}\widehat{\Phi}(\breve x^k,\breve y^k) \end{pmatrix} \] are formed. For \( \eta_k:=s_k^\top q_k>0 \), one sets \( \alpha_{\grave{k}} = \tilde{\eta}_k\, s_k^\top s_k \), and subsequently restricts it to an interval \( [\alpha_{\min},\alpha_{\max}] \). If \(\eta_k\le 0\), the previous usable step size is retained, or a conservative initial step size is used.

Algebraic CO

Purpose.
The Merit phase is a first-order method. It reduces residuals and the duality gap, but does not always deliver the active basis exactly. Therefore, if required, an algebraic CO follows. It attempts to reconstruct an exactly feasible primal-dual solution from the approximation point \((x,y)\).

Active candidates.
Primal active constraints have those indices \(i\) for which \(a_i^\top x\approx b_i\) holds. Active non-negativity conditions have those indices \(j\) for which \(x_j\approx0\) holds. From these candidates, linear systems of size \(n\) are formed. A primal candidate is accepted if the reconstructed \(x\) satisfies \[ x\ge0,\qquad Ax\le b. \] A dual candidate \(y\) is then sought with \[ y\ge0,\qquad A^\top y\ge c. \] Only if, in addition, \(|b^\top y-c^\top x|\) is sufficiently small is the pair accepted as an optimal solution.

Certified CO.
An isolated good primal candidate is not sufficient. Likewise, an isolated good dual candidate is not sufficient. Only a primal-dual pair satisfying simultaneously \[ p_\infty=0,\qquad d_\infty=0,\qquad g_{\mathrm{abs}}=0 \] up to the required tolerance is accepted. This prevents a merely formally feasible dual candidate with an incorrect objective value from replacing the solution.

Staged CO.
The staged CO substantially limits the combinatorial cost.

In Stage 1, the probably active constraints are obtained from small slacks \(b_i-a_i^\top x\) and variables \(x_j.\) This stage is fast, but may fail if the active sets are poorly separated.

In Stage 2, the candidate set is enlarged. Additional constraints with medium-sized slacks are admitted. This increases the probability of success, but also the number of tested bases.

In Stage 3, a complete or nearly complete fallback is admitted. This stage is expensive, but robust. It is used only if the previous stages do not produce a certified primal-dual solution.

The CO can also be switched off completely. In that case, the method outputs the best Intex–Merit point found so far, without algebraic basis reconstruction. This algorithmic limitation is accompanied by a strictly memory-efficient implementation at the software level. The inner loop of the combinatorial candidate search operates entirely without allocation.

All vectors and index arrays are pre-allocated, and only overwritten during each iteration. This prevents any dynamic memory pressure and ensures that the garbage collection of the runtime environment is not triggered during the critical phase. Due to this mechanical efficiency, the theoretical worst-case complexity of the basis search loses its impact in practice, as traversing the permutations only requires minimal CPU cycles.

Timeout.
If a time limit is reached, the method terminates in a controlled way and returns the best state known up to that point: \[ x_{\mathrm{best}},\qquad y_{\mathrm{best}},\qquad p_\infty,\qquad d_\infty,\qquad g_{\mathrm{abs}},\qquad p^*,\qquad d^*. \] If no certified solution has been found up to that point, this is explicitly indicated. A timeout is therefore not a mathematical false conclusion, but an incomplete computational state.

A primal infeasible LP has no \(x\ge0\) with \(Ax\le b\): neither can Intex lower the radius to \(r=0\), nor can Merit enforce \(p_\infty=0\). Typically, \( p_\infty>0\) remains, or the radius reduction stagnates. If the dual LP is infeasible, then there is no \(y\ge0\) with \(A^\top y\ge c.\) Usually \(d_\infty>0.\) In primal terms, this often corresponds to an unbounded primal problem.

If the problem is formally solvable but numerically badly scaled or degenerate, Intex and Merit may stagnate even though a solution exists. The CO then decides whether the active structure can nevertheless be reconstructed algebraically. If this fails, the method does not provide a false proof of optimality, but a non-certified approximation point.

Thus the method distinguishes, in practice, four cases: \[ \begin{array}{ll} p_\infty=d_\infty=g_{\mathrm{abs}}=0 &\text{certified as solved},\\[0.3em] p_\infty>0 &\text{primal feasibility not reached},\\[0.3em] d_\infty>0 &\text{dual feasibility not reached},\\[0.3em] p_\infty,d_\infty\approx0,\ g_{\mathrm{abs}}>0 &\text{feasible, but optimality not reached}. \end{array} \] A complete infeasibility certificate in the sense of a Farkas witness is therefore not automatically included. However, it may be added as an additional diagnostic stage.

Structured linear algebra by FTD and FSR

Motivation.
Intex, Merit, and CO repeatedly require operations of the form \[ Ax,\qquad A^\top y,\qquad b^\top y,\qquad c^\top x. \] For a general sparse matrix this costs \( \mathcal O(\operatorname{nnz}(A)).\) For a dense unstructured matrix, the cost is \(\mathcal O(mn).\) If \(A\) has structure, this need not be the case. For Toeplitz, Hankel, block-circulant, or convolution-like matrices, the operator evaluations can be replaced by FTD and FSR.

FTD.
In practice, FTD can often describe matrices with \(mn\) entries by only a few generators, diagonal sequences, displacement defects or convolution data, and can therefore treat them operatorially. A suitable model form is\[ A\approx K_s+UV^\top+E, \qquad K_s=T_1T_2\cdots T_s. \] Here \(T_1,\dots,T_s\) are Toeplitz, Hankel, block-circulant or related convolution operators, \[ U,V\in{}^{\nu}\mathbb R^{n\times r}, \qquad r\ll n, \] describe a low-rank boundary or defect correction, and \(E\) denotes the remaining approximation error. The case \(s=1\) yields a simple Toeplitz or convolution kernel. The case \(s>1\) permits a multiply iterated generation of the structured kernel.

The advantage of this multiple iteration is that the remaining numerical rank of the correction may decrease. Formally, this can be expressed by \[ r_\varepsilon(s) = \min_{T_1,\dots,T_s} \operatorname{rank}_\varepsilon \left( A-T_1T_2\cdots T_s \right) \] Here \(r_\varepsilon(s)\) denotes the numerical residual rank after the choice of \(s\) structured factors. An additional factor is worthwhile in particular when the reduction of \(r_\varepsilon(s)\) outweighs the additional cost of its evaluation.

The factors are obtained from local averages, diagonal sums, displacement defects and recursive balancing steps. Instead of the full matrix, only generator, diagonal, convolution and correction data are stored. Matrix-vector products, transposed products and selected projections can thereby be carried out structurally.

If the effective structural width is \(\ell\) and \( p=\ell n, \) then the structured kernel computation typically has cost of the form \( \mathcal O(p^2) \) rather than dense unstructured cost. For FFT-capable convolution kernels, matrix-vector products can also be organised in the form \( \mathcal O(s\,n {}_2n+nr) \) where \(s\) denotes the number of structured factors and \(r\) the rank of the correction.

FSR.
The FSR complements the FTD arithmetically. While the FTD provides a structured representation of the operator, FSR organises the actual computation recursively and in a word-length-preserving manner. If \(\ll v\) denotes a shift by \(v\) digits, a matrix or an operator can be decomposed into high- and low-weight parts: \[ A=(A_{\mathrm{hi}}\ll v)+A_{\mathrm{lo}}, \] Products, Schur complements, inverse subsystems and reconstructions are built recursively from smaller parts. Thus intermediate values remain controlled, and BigInt or BigFloat arithmetic does not grow unnecessarily. The recursion is well suited to parallelisation here. Hence FSR is particularly suitable for high precision, exact integer arithmetic and large structured linear systems.

Low-rank corrections.

If \( A\approx K_s+UV^\top \) and the kernel solves with \(K_s\) are structurally available, the low-rank correction can additionally be split off by the Woodbury formula: \[ (K_s+UV^\top)^{-1} = K_s^{-1} – K_s^{-1}U \left(I+V^\top K_s^{-1}U\right)^{-1} V^\top K_s^{-1}. \] This reduces a large linear system to structured kernel solves and a small correction system. For the simple Toeplitz kernel \( A=T+UV^\top \) one furthermore obtains the resolvent form \[ \det(\lambda I-A) = \det(\lambda I-T) \det\!\left(I_r-V^\top(\lambda I-T)^{-1}U\right). \] The second eigenproduct has size only \(r\times r\). Eigenvalue components caused by the boundary or defect correction can therefore be investigated as zeros of the small rational equation \[ \det\!\left(I_r-V^\top(\lambda I-T)^{-1}U\right)=0 \] Thus outlier eigenvalues can be separated from the kernel spectrum determined by \(T\).

Interaction.
The FTD provides the structured decomposition \(A\approx T_1T_2\), while the FSR performs the resulting products, inversions, Schur-complement steps, and back-transformations recursively. Here the combination replaces the general operator evaluation \(\mathcal O(dmn)\) by \(\mathcal O(p^2).\)

Integration into Intex–Merit.
The FTD/FSR structure accelerates the computation in three ways:

First, the Intex phase, because bounds, cuts and centroid movements are repeatedly computed there by operations with \(A\) and \(A^\top\).

Second, the Merit polish, because the function value and gradient essentially consist of: \[ Ax-b,\qquad c-A^\top y,\qquad b^\top y-c^\top x. \] Third, possibly the CO. Solving a dense active system of size \(n\) without structure costs \( \mathcal O(n^3). \) If, however, the active subsystem has the same Toeplitz, Hankel, block-circulant or convolution-like structure, FTD/FSR can reduce this step to structured kernel operations, Schur-complement steps and small correction systems.

Thus Intex, Merit, CO, FTD, and FSR do not form separate building blocks, but rather a coherent architecture: Intex = geometric reduction,\(\quad\) Merit = numerical polishing,\(\quad\) CO = algebraic certification\(\quad\) FTD/FSR = structured linear algebra.

Let \(\mathcal M_A\) denote the cost of a complete structure-preserving operator evaluation with \(A\) and \(A^\top\). Depending on the representation of \(A\), one obtains \[ \mathcal M_A= \begin{cases} \mathcal O\!\left(\operatorname{nnz}(A)\right), & \text{for sparse matrices},\\[0.3em] \mathcal O(mn), & \text{for dense unstructured matrices},\\[0.3em] \mathcal O(p^2), & \text{for }\mathrm{FTD}/\mathrm{FSR}\text{-structured kernel operations},\\[0.3em] \mathcal O\!\left(s\,n\,{}_2n+nr\right), & \text{for }\mathrm{FFT}\text{-capable convolution factors}\\ & \text{with rank-}r\text{ correction}. \end{cases} \] Here \(p=\ell n\) denotes the effective structured problem size for structural width \(\ell\), \(s\) the number of structured factors in the kernel \( K_s=T_1T_2\cdots T_s, \) and \(r\) the rank of the low-rank boundary or defect correction. The Intex phase therefore has order \( \mathcal O\!\left({}_2(\tilde\alpha\rho)^2\,\mathcal M_A\right). \) For \(K\) accepted iterations and \(B\) additional backtracking evaluations, the Merit polish requires the cost \( \mathcal O\!\left((K+B)\mathcal M_A\right). \)

For the CO without structure, without additional structure, \( \mathcal C=\mathcal O(Tn^3), \) if \(T\) active basis candidates are actually tested. Since, for \(Q\) candidate equations, the unfavourable case may satisfy \( T\le \binom{Q}{n} \) an uncontrolled CO is potentially combinatorial. The staged model therefore first reduces \(Q\) and activates larger candidate sets only when smaller active systems do not provide sufficient certification.

If the active subsystem has the same Toeplitz, Hankel, block-circulant or convolution-like structure as the original operator, the CO can be carried out by structured kernel operations, Schur-complement steps and small correction systems rather than by dense factorisations. Then approximately \( \mathcal C=\mathcal O(Tp^2) \) holds or, for FFT-capable product kernels and low-rank correction, a corresponding dependence on \( s\,n{}_2n+nr. \)

The total order of the hybrid method is therefore \[ \mathcal O\!\left( {}_2(\tilde\alpha\rho)^2\mathcal M_A + (K+B)\mathcal M_A + \mathcal C \right). \] The actual runtime is thus determined by three quantities: the number of geometric Intex steps, the number of Merit evaluations and the cost of algebraic certification in the CO

In addition to the matrix or generator representation, the memory requirement in the Intex–Merit kernel is \( \mathcal O(m+n). \) With dense storage this gives altogether \( \mathcal O(mn+m+n)=\mathcal O(mn). \) With structured storage, the memory requirement depends on the generator width, the number of factors and the rank of the correction. For a representation \( A\approx K_s+UV^\top+E \) it is closer to \( \mathcal O(p+nr+m+n), \) where \(p\) describes the structured kernel representation and \(nr\) the storage of the low-rank correction.

Corollary from the INLIF method.

Every LP can be solved in \(\mathcal O(Ln^3)\) or certified as infeasible.\(\square\)

The LINEAL method (Linear Interpolation, Nearing Extrapolation, and Amending Lifting) combines Intex and INLIF into a unified solution approach for linear inequality systems and linear programmes. Intex approaches a solution through linear interpolation and extrapolation within successively reduced relaxations; any remaining or stagnating violations may then be absorbed and back-substituted by the amending lifting mechanism of INLIF. Intex therefore remains an independently effective primary method, while INLIF serves, where required, as an algebraic completion procedure.

A second solution on the same objective face.
If \(x^o\) is an optimal solution, the LP \[ \max \left\{ \|x-x^o\|_1: c^\top x=c^\top x^o,\; Ax\le b,\; x-x^o\in[-1,1]^n,\; x\in{}^{\nu}\mathbb R_{\ge0}^n \right\} \] can determine a second optimal solution, provided that the optimal face contains more than one point. The dual vector \(y^o\) can be treated analogously.

Linear systems.
The equivalent representation \(Ax\le b,\ -Ax\le -b\) turns a linear system of equations \(Ax=b\) into an INLIF system, which can be treated by means of the merit function \(\widehat{\Phi}(x)=\|Ax-b\|_2^2\). In the regular case this leads to the solution of the system; in the singular or inconsistent case, only to a least-squares solution. The gradient is \(\nabla\Phi(x)=A^\top(Ax-b)\). Additional side constraints such as \(x_j\ge0\) can be incorporated by projection or by an LP formulation.

Thus Gaussian elimination is not replaced, but embedded into a feasibility-oriented reduction view. The advantage may lie in early termination: INLIF can stop as soon as the transformed origin becomes feasible or a residual-cone row provides an infeasibility certificate. For regular systems of equations, Gaussian elimination remains the natural direct method; for overdetermined, inconsistent, or inequality-coupled systems, however, INLIF is plausible as a structural alternative.

The reduction logic of INLIF is violation-driven. The current tableau distinguishes between violated activation rows and satisfied passive rows. Satisfied rows remain part of the system, but do not force an absorption step. This allows the method to restrict the current reduction to the violated constraints and to defer non-violated inequalities until their possible later activation.

Regularity.
The matrix \(A\) is regular if and only if the homogeneous system \(Ax=0\) has only the zero solution. This can be tested by an LP of the form \[ \max\{\|x\|_1:Ax=0,\;x\in[-1,1]^n\} \] If the optimal value is zero, the kernel is trivial.

Inverse columns.
For regular \(A\in{}^{\nu}\mathbb R^{n\times n}\), each column \(\alpha_j\) of \(A^{-1}\) can be determined by \(A\alpha_j=e_j\), where \(e_j\) is the \(j\)-th unit vector. Without structure, this is a classical linear system; with FTD/FSR, the structure of \(A\) can be exploited in the solution.

Eigenvalue problems.
Eigenvalue equations \(Ax=\lambda x\) can be transferred into constraints of an extended optimisation problem. In practice, this is of interest when additional constraints, normalisations, or structural conditions are imposed on \(x\) and \(\lambda\). For general dense matrices, classical spectral methods usually remain more suitable; however, for structured or constrained formulations, the Merit-LP viewpoint can be advantageous.

Polynomial and Padé approximation.
Least-squares problems of the form \[ y=r+Xc,\qquad X^\top r=0 \] can also be formulated as linear or quadratic Merit problems. Additional constraints on the coefficients \(c\) can be inserted directly. For Padé approximations, the discrete Fourier transform can provide structured convolution components, which are then processed by FTD/FSR.

Convex programmes.
For convex programmes \(\quad \min \left\{ f_1(x): x\in{}^{\nu}\mathbb R^n,\; (f_2(x),\dots,f_m(x))^\top\le0 \right\}\quad\) the same idea can be transferred to nonlinear residuals. The linear terms \(Ax-b\) are then replaced by the violations \(f_k(x)\). The basic structure remains: \[ \text{feasibility}+\text{duality or optimality condition} \quad\longrightarrow\quad \text{Merit function}. \] A polynomial runtime statement, however, additionally depends on the evaluability, smoothness, conditioning, and structure of the functions \(f_k\).

Shift-polynomial FSR multiplication for integer matrices

For large entries, the classical multiplication of two matrices over \(\mathbb Z\) is determined not only by the number of scalar products, but also substantially by the cost of integer multiplication. For \( A,B\in \mathbb Z^{N\times N} \) with entries of large bit length, a recursive block structure may therefore be advantageous if it improves memory locality, arranges the partial products in a favourable accumulation order, and avoids temporary product matrices.

In what follows, assume first that \(N\) is divisible by \(3\). Write \[ A= \begin{pmatrix} A_{11}&A_{12}&A_{13}\\ A_{21}&A_{22}&A_{23}\\ A_{31}&A_{32}&A_{33} \end{pmatrix}, \qquad B= \begin{pmatrix} B_{11}&B_{12}&B_{13}\\ B_{21}&B_{22}&B_{23}\\ B_{31}&B_{32}&B_{33} \end{pmatrix}, \] where all blocks have the same size.

The shift-polynomial approach formally encodes the three block columns of \(A\) and the three block rows of \(B\) by powers of a central shift symbol \(s\). For the three inner couplings one defines \[ P_1 = (A_{11}+sA_{21}+s^2A_{31}) (B_{11}+s^3B_{12}+s^6B_{13}), \] \[ P_2 = (A_{12}+sA_{22}+s^2A_{32}) (B_{21}+s^3B_{22}+s^6B_{23}), \] \[ P_3 = (A_{13}+sA_{23}+s^2A_{33}) (B_{31}+s^3B_{32}+s^6B_{33}). \] Then \[ P=P_1+P_2+P_3 \] is a formal matrix polynomial in \(s\), whose coefficients directly contain the resulting blocks of \(AB\). More precisely, \[ (AB)_{pq} = [s^{\,p-1+3(q-1)}]P \qquad (p,q\in\{1,2,3\}). \] Thus the order of the coefficients is \[ s^0,s^1,s^2,s^3,s^4,s^5,s^6,s^7,s^8 \] and corresponds to the block positions \[ 11,21,31,12,22,32,13,23,33. \]

Correctness theorem

Let \(A,B\in R^{N\times N}\) be matrices over an associative, in general non-commutative ring \(R\). Let the shift symbol \(s\) be central, i.e. \( sX=Xs \) for all occurring blocks \(X\). Then the coefficient extraction above yields exactly the matrix product \(AB\).

\emph{Proof.} For \(r\in\{1,2,3\}\), the expression \[ (A_{1r}+sA_{2r}+s^2A_{3r}) (B_{r1}+s^3B_{r2}+s^6B_{r3}) \] contains the term \( A_{pr}B_{rq} \) exactly at the exponent \( (p-1)+3(q-1). \) After summation over \(r=1,2,3\), the coefficient of \( s^{p-1+3(q-1)} \) therefore gives \( {\LARGE{\textbf{+}}}_{r=1}^3 A_{pr}B_{rq} = (AB)_{pq}. \) The order of the factors is preserved; in particular, no commutativity of the matrix blocks is assumed.\(\square\)

Recursive FSR form

If the construction above is applied recursively to the block products, one obtains a ternary FSR multiplication. The recursion need not be carried down to the scalar level. For the leaf size \(L\), a hybrid variant is more favourable: \[ N \longrightarrow \frac{N}{3} \longrightarrow \frac{N}{3^2} \longrightarrow \cdots \longrightarrow L. \] At the leaf level, the method returns to classical integer matrix multiplication.

The recursion depth is therefore \( d={}_{3}(N/L), \) provided that \(N/L\) is a power of \(3\). For \( N=729,\) and \(L=27 \) one obtains, for example, \( d=3. \) However, the method is not restricted to pure powers of three. It is sufficient that \( N=3^dL \) with a practically suitable leaf size \(L\). For example, \( 1026=3^3\cdot38,\) and \( 1053=3^3\cdot39. \) The choice of the leaf size is essential. A leaf size that is too small leads to deeper recursion, many small leaf calls, and increased administrative and memory overhead.

The practical operation is not \( \qquad T=AB,\qquad C\leftarrow C+T,\qquad \) but directly \( C\leftarrow C+AB. \) This avoids temporary product matrices. For large \(\mathbb Z\)-entries this is essential, since memory management can otherwise account for a substantial part of the running time.

Signs

Since the recursive form accumulates directly into the target matrix and does not recover coefficient windows by integer divisions, the method is also correct for mixed-sign BigInt entries. The signs are handled entirely by the ordinary multiplication and addition in \(\mathbb Z\). This distinguishes the present form from pure packing methods, in which several coefficients are encoded in one large integer and subsequently separated again by shifts or divisions. For mixed signs, such recovery requires additional balancing or division steps. The FSR form used here avoids this recombination.

Empirical findings

The following measurements are empirical running-time results for dense BigInt matrices. They depend on processor architecture, thread count, Julia version, GMP binding, leaf size, and memory management. Nevertheless, they clearly show that the recursive administrative overhead can be more than compensated for by sufficiently large matrices.

For dense matrices of size \( 2187\times2187 \) with non-negative BigInt entries of \(128\) bits and hybrid recursion down to the leaf size \( L=81 \) a serial speedup of approximately \( 3.54 \) over classical GMP-based matrix multiplication was obtained on a \(24\)-thread workstation. In the parallel variant, a speedup of approximately \( 1.38 \) over the classically parallelised multiplication was measured.

For dense matrices of size \( 729\times729 \) with non-negative BigInt entries of \(256\) bits and hybrid recursion down to the leaf size \( L=27 \) a serial speedup of approximately \( 3.36 \) over classical GMP-based matrix multiplication was obtained on the same workstation. In the parallel variant, a speedup of approximately \( 1.33 \) over the classically parallelised multiplication was measured.

The effect also persists for mixed-sign entries. For dense matrices of size \( 729\times729 \) with mixed \(256\)-bit BigInt entries and leaf size \( L=27 \) a serial speedup of approximately \( 3.80 \) and a parallel speedup of approximately \( 1.41 \) over the corresponding classical GMP-based multiplication were measured on an \(24\)-thread system.

For the size \( 1026\times1026 \) with mixed \(256\)-bit BigInt entries one has \( 1026=3^3\cdot38. \) With leaf size \( L=38 \) a serial and parallel speedup of approximately \( 4.22 \) and \( 1.32 \) were obtained. This shows that the method is not restricted to pure powers of three, but remains useful for sizes of the form \(3^dL\). The optimal leaf size is not necessarily the largest still reasonable leaf size. Smaller leaves increase recursion and memory overhead, but can nevertheless be faster because of more favourable leaf multiplications.

Comparison with binary FastShift

A neighbouring practical comparison with a binary FastShift method yielded, for mixed \(256\)-bit BigInt entries of size \( 1024\times1024 \) on \(24\) threads, a running time of approximately \( 24.84\ \mathrm{s}. \) The ternary shift-polynomial FSR form, by contrast, computed the nearby product of size \( 1026\times1026 \) in approximately \( 5.30\ \mathrm{s}. \) Although the matrix sizes are not identical, this comparison indicates a practical advantage of the ternary form in the setting considered. An essential reason is that the binary packing method works with wider intermediate numbers and subsequent recombination, whereas the ternary FSR form accumulates directly into the target blocks.

Classification

Shift-polynomial FSR multiplication is not to be understood as a replacement for highly optimised floating-point BLAS multiplication. Its natural area of application is rather exact integer arithmetic, rational arithmetic, and structured recursions in which the costs of the individual arithmetic operations are high.

The advantage does not arise from a mere reformulation of matrix multiplication, but from the combination of \[ \text{large matrix size} + \text{large integer arithmetic} + \text{hybrid ternary recursion} + \text{direct in-place accumulation}. \] The choice of the leaf size is particularly important. A recursion that is too deep, for instance with a very small leaf size, increases the administrative and memory overhead. A recursion that is too shallow, on the other hand, does not make sufficient use of the structure. The optimal leaf size is therefore an algorithmic parameter and depends on the matrix size, bit length, processor architecture, and parallelisation.

Three-branch projector structure for linear systems

Basic block.

For \(a,c,d\in\mathbb K\), \(\mathbb K\in\{\mathbb R,\mathbb C\}\), define \[ M(a,c,d) = \begin{pmatrix} \check{(a+d)}&0&\check{(a-d)}\\ 0&c&0\\ \check{(a-d)}&0&\check{(a+d)} \end{pmatrix}. \] For a vector \(x=(x_1,x_2,x_3)^\top\), introduce the local branch coordinates \[ s=\check{(x_1+x_3)},\qquad m=x_2,\qquad q=\check{(x_1-x_3)}. \] The inverse transformation is \[ x_1=s+q,\qquad x_2=m,\qquad x_3=s-q. \] In these three branches, \(M(a,c,d)\) acts diagonally: \[ s\mapsto as,\qquad m\mapsto cm,\qquad q\mapsto dq. \] In particular, \[ M(a,c,d)M(e,f,g)=M(ae,cf,dg), \] and \[ M(1,1,1)=I_3. \] The class is therefore closed under addition and multiplication and forms a commutative algebra isomorphic to \(\mathbb K^3\). For \(a,c,d\in\mathbb K^*\), \[ M(a,c,d)^{-1}=M(\tilde a,\tilde c,\tilde d). \]

Block systems of size \(3m\).
Let \[ A^{(s)},A^{(m)},A^{(q)}\in\mathbb K^{m\times m}. \] From these three branch matrices, construct a matrix \[ A_0\in\mathbb K^{3m\times3m} \] whose \(3\times3\) block at position \((i,j)\) is given by \[ (A_0)_{ij} = M\!\left(a^{(s)}_{ij},a^{(m)}_{ij},a^{(q)}_{ij}\right). \] Thus \(A_0\) has a three-branch projector structure.

For a linear system \( A_0x=b, \) decompose \(x\) and \(b\) blockwise into the branches \[ s_j=\check{(x_{3j-2}+x_{3j})},\qquad m_j=x_{3j-1},\qquad q_j=\check{(x_{3j-2}-x_{3j})}, \] and correspondingly \[ b^{(s)}_j=\check{(b_{3j-2}+b_{3j})},\qquad b^{(m)}_j=b_{3j-1},\qquad b^{(q)}_j=\check{(b_{3j-2}-b_{3j})}. \] Then the original \(3m\times3m\) system splits exactly into three ordinary \(m\times m\) systems: \[ A^{(s)}s=b^{(s)},\qquad A^{(m)}m=b^{(m)},\qquad A^{(q)}q=b^{(q)}. \] After solving these systems, \(x\) is reconstructed by \[ x_{3j-2}=s_j+q_j,\qquad x_{3j-1}=m_j,\qquad x_{3j}=s_j-q_j. \] Consequently, \[ A_0\text{ is regular} \quad\Longleftrightarrow\quad A^{(s)},A^{(m)},A^{(q)}\text{ are regular}. \]

Schur extension.
For system sizes that are not fully described by the three-branch structure, let \[ A= \begin{pmatrix} A_0&R\\ R^\top&C \end{pmatrix}, \] where \( A_0\in\mathbb K^{3m\times3m} \) has the three-branch projector structure, while \( R\in\mathbb K^{3m\times r} \) and \( C\in\mathbb K^{r\times r} \) describe the residual and boundary couplings. The linear system \[ \begin{pmatrix} A_0&R\\ R^\top&C \end{pmatrix} \begin{pmatrix} x\\y \end{pmatrix} = \begin{pmatrix} b_0\\b_1 \end{pmatrix} \] is solved by means of the Schur complement. From \( A_0x+Ry=b_0 \) one obtains \( x=A_0^{-1}(b_0-Ry). \) Substitution into \[ R^\top x+Cy=b_1 \] gives \[ \left(C-R^\top A_0^{-1}R\right)y = b_1-R^\top A_0^{-1}b_0. \] The Schur complement is therefore \[ W=C-R^\top A_0^{-1}R. \] After solving the small system \[ Wy=b_1-R^\top A_0^{-1}b_0, \] the structured part is reconstructed by \[ x=A_0^{-1}(b_0-Ry). \] Here \(A_0^{-1}\) is not formed explicitly. The products \[ A_0^{-1}R,\qquad A_0^{-1}b_0 \] are computed through the three branch systems \[ A^{(s)},\qquad A^{(m)},\qquad A^{(q)}. \]

SPD case.
If \[ A^{(s)},A^{(m)},A^{(q)} \] are symmetric positive definite and the Schur complement \[ W=C-R^\top A_0^{-1}R \] is also symmetric positive definite, then the entire computation can be carried out by Cholesky factorisations: \[ A^{(s)}=L_sL_s^\top,\qquad A^{(m)}=L_mL_m^\top,\qquad A^{(q)}=L_qL_q^\top, \] and \[ W=L_WL_W^\top. \] This is numerically preferable to a general dense factorisation.

Complexity.
A direct dense solution of a system of size \( 3m+r \) has cubic order \( \mathcal O\!\left((3m+r)^3\right). \) In the SPD case, the structured Schur variant requires three factorisations of size \(m\times m\), the construction of the Schur complement and the solution of the small \(r\times r\) system. This gives approximately \[ \mathcal O\!\left( 3m^3 + 3m^2r + mr^2 + r^3 \right). \] For fixed small \(r\), the dominant term is \( \mathcal O(3m^3) \) instead of \( \mathcal O(27m^3). \) The ideal cubic structural gain is therefore approximately \[ \frac{(3m)^3}{3m^3}=9. \] The coupling terms \( A_0^{-1}R \) and \( R^\top A_0^{-1}R \) reduce this ideal value in practice, but do not remove the main asymptotic advantage for \(r\ll m\).

Numerical example.
In a Julia test, a symmetric positive definite system of size \( 3m+r \) was generated. The structured core \(A_0\) consisted of \(3\times3\) blocks of the form \(M(a,c,d)\), while \(R\) described a small boundary coupling and \(C\) a small positive definite residual block.

For \( m=100,\qquad r=5, \) one obtains \( 3m+r=305. \) The direct Cholesky solution required approximately \( 0.001232\,\mathrm{s}, \) whereas the structured Schur–Cholesky solution required approximately \( 0.000285\,\mathrm{s}. \) The errors were of the order \[ \|x_{\mathrm{Schur}}-x_{\mathrm{true}}\|_2 \approx 7.909\cdot10^{-15}, \] and \[ \|Ax_{\mathrm{Schur}}-b\|_2 \approx 5.414\cdot10^{-14}. \]

For \( m=300,\qquad r=5, \) one obtains \( 3m+r=905. \) The direct Cholesky solution required approximately \( 0.006438\,\mathrm{s}, \) whereas the structured Schur–Cholesky solution required approximately \( 0.004157\,\mathrm{s}. \)

For \( m=500,\qquad r=5, \) one obtains \( 3m+r=1505. \) The direct Cholesky solution required approximately \( 0.052780\,\mathrm{s}, \) whereas the structured Schur–Cholesky solution required approximately \( 0.008407\,\mathrm{s}. \) Thus, in the largest test, a practical speed-up factor of approximately \[ \frac{0.052780}{0.008407}\approx 6.28 \] was obtained. The theoretical cubic comparison gives, in this case, \[ \frac{1505^3}{3\cdot500^3+5^3}\approx9.09. \] The deviation from the ideal value is explained by coupling costs, memory access and implementation constants.

Interpretation.
The method does not accelerate arbitrary dense linear systems. Its advantage appears precisely when the leading core has a three-branch projector structure, or when it approximately assumes such a structure after a suitable permutation or transformation. In this case, a large system is not treated as one monolithic dense linear system, but is split into several smaller branch systems and a small Schur complement: \[ \boxed{ 3m+r \quad\longrightarrow\quad 3\cdot m+r. } \] If the branch systems are recursively nested, this decomposition can be continued. For \(k\) three-branch levels, one formally obtains a decomposition into \( 3^k \) independent branches. Hence the method fits naturally with recursive structured procedures such as FSR.

Let

be a polynomial map of degree at most \(2\). Suppose that

for all \((x,y)\in\mathbb{C}^2\). The theorem treats the quadratic special case of the Jacobian conjecture. It is shown that \(H\), being a polynomial automorphism, has no genuine double-point pair.

Lemma on the separation of infinitesimal orders

Let \(\delta:=\tilde{\omega}\ne0\) be infinitesimal, and let \(a_r,\ldots,a_m\) be finite quantities independent of \(\delta\). If

then

Proof.

Suppose that not all coefficients vanish. Let \(j\) be the smallest index for which \(a_j\ne0\). Equation (1) then yields

Since \(\delta\ne0\), it follows that

The right-hand side is infinitesimal. By contrast, the coefficient \(a_j\), which is independent of \(\delta\), is either zero or non-infinitesimal. Since \(a_j\ne0\), this is a contradiction.

Hence all coefficients vanish.\(\square\)

Meaning of the order lemma

In particular, the lemma states that terms of different infinitesimal orders cannot cancel one another. A term of order \(\delta^3\) cannot cancel a non-zero term of order \(\delta^2\). More generally,

is possible only if \(a_r=0\) first. The same argument is then applied to the next higher order.

This is the infinitesimal form of comparing homogeneous polynomial degrees.

Degree zero

If \(H\) has degree zero, then \(H\) is constant. Hence \( JH=0 \) and therefore \( \operatorname{eig} JH=0. \) This contradicts \(e\ne0\). Consequently, no Keller map of degree zero exists.\(\square\)

Degree one

If \(H\) has degree one, then \(H\) has the form \( H(z)=c+Az \) with \( z=(x,y)^\top, c\in\mathbb{C}^2 \) and \(A\in\mathbb{C}^{2\times2}\).

Since \( JH=A \) and \( \operatorname{eig} A=e\ne0, \) the matrix \(A\) is invertible, with inverse map

Thus \(H\) is an affine automorphism and has no genuine double-point pair.\(\square\)

Normalisation of the quadratic case

Now let \(H\) have degree exactly \(2\). Then \(H\) has a unique decomposition

where \(Q\) is a homogeneous quadratic vector polynomial.

At the origin, \( JH(0)=A. \) The invertibility of \(A\) therefore follows from the Keller condition \( \operatorname{eig} A=e\ne0. \) Define

Then, with the homogeneous quadratic polynomial \( R(z):=A^{-1}Q(z), \) one has

Since \( JK(z)=A^{-1}JH(z), \) it follows that

Separation of the linear and quadratic eigenproduct parts

Since \(R\) is homogeneous quadratic, \( M(z):=JR(z) \) is a homogeneous linear matrix function. One has \( JK(z)=I+M(z). \) For a \(2\times2\) matrix \(M\),

The Keller condition therefore gives

Now replace \(z\) by \(\delta z\). Homogeneity gives \( M(\delta z)=\delta M(z), \) and hence

The first summand is of order \(\delta\), whereas the second is of order \(\delta^2\). By the order lemma, they cannot cancel one another. Division by \(\delta\ne0\) yields successively \( \operatorname{tr}M(z)=0 \) and \( \operatorname{eig} M(z)=0. \) Since \(z\) was arbitrary, the polynomial identities \( \operatorname{tr}JR=0 \) and \( \operatorname{eig} JR=0 \) hold. This explicitly takes into account that the quadratic infinitesimal term cannot cancel the linear one.

Structure of the quadratic part

Write \( R=(P,S), \) where \(P\) and \(S\) are homogeneous quadratic polynomials.

From

it follows that there exists a homogeneous cubic polynomial \(h\) such that

Thus

Together with \( \operatorname{eig} JR=0, \) this gives

Write

Then

Substitution into equation (12), followed by comparison of coefficients, yields

If \(a\ne0\), set \( r:=\tilde a b. \) Then \( c=ar^2 \) and \( d=ar^3. \) Consequently,

If \(a=0\), then \(ac=b^2\) first gives \(b=0\), and \(bd=c^2\) subsequently gives \(c=0\). Hence

In both cases, there exist a linear form \( \ell(x,y)=\alpha x+\beta y \) and a constant \(\kappa\in\mathbb{C}\) such that \( h=\kappa\ell^3. \) Equation (11) then yields

With \( w:=3\kappa(\beta,-\alpha), \) one has \( R(z)=w\ell(z)^2. \) Moreover,

Thus every normalised quadratic Keller map has the form

Direct exclusion of double points

Suppose that there exist \(z_1\ne z_2\) such that

Applying the linear form \(\ell\) to \(K(z)\), and using \(\ell(w)=0\), gives

Equation (19) therefore implies

Hence the quadratic correction terms coincide:

It then follows from

that \( z_1=z_2. \) This contradicts the assumption \(z_1\ne z_2\). Consequently, \(K\) has no genuine double-point pair.\(\square\)

Explicit inverse map

The normalised map \( K(z)=z+w\ell(z)^2 \) has the polynomial inverse

Indeed, \(\ell(w)=0\) gives

Therefore,

For the original map \( H(z)=c+A K(z), \) set \( u:=A^{-1}(y-c). \) Then the inverse map is polynomial and is given by

Canonical normal form

If \(w=0\), then \(K\) is the identity. Hence assume \(w\ne0\). Choose a vector \(r\in\mathbb{C}^2\) such that \( \ell(r)=1. \) Since \(\ell(w)=0\), the vectors \(w\) and \(r\) form a basis of \(\mathbb{C}^2\). Writing \( z:=uw+vr, \) one has \( \ell(z)=v \) and hence

In the coordinates \((u,v)\), the map is

Up to affine coordinate changes, every quadratic Keller map therefore has one of the two forms

The Keller class \(F=p(x+f(x,y))\) and \(G=q(y+f(x,y))\)

Let \( f\in\mathbb C[x,y] \) and \( p,q\in\mathbb C[t]. \) Consider the polynomial map

with

and

Suppose that the Keller condition

holds for all \((x,y)\in\mathbb C^2\), where \( e\in\mathbb C^*. \)

It will be shown that the Keller condition necessarily reduces the outer polynomials \(p\) and \(q\) to affine polynomials and that \(f\) can be determined completely. An explicit polynomial inverse map will then be given.

Introduction of the inner coordinates

Set

and

The partial derivatives of \(u\) and \(v\) are

and

Calculation of the Jacobi matrix

By the chain rule,

and

Thus the Jacobi matrix has the form

Its eigenproduct is

Hence the fundamental eigenproduct identity is

The Keller condition is therefore equivalent to

Unit argument in the polynomial ring

The polynomial ring \( \mathbb C[x,y] \) is an integral domain. Its units are precisely the non-zero constant polynomials. The right-hand side of equation (37) is a unit in \(\mathbb C[x,y]\). Consequently, each of the three factors on the left-hand side must itself be a unit. Hence there exist constants

such that

and

Necessarily,

Non-constancy of the inner polynomials

The polynomial \(u=x+f(x,y)\) cannot be constant.

Suppose that \( u=\gamma \) for some constant \(\gamma\in\mathbb C\). Then \( F=p(\gamma) \) would be constant. Consequently, the first row of the Jacobi matrix of \(H\) would vanish, and one would have \( \operatorname{eig} JH=0. \) This contradicts the assumption \(e\ne0\). Thus \( u=x+f(x,y) \) is non-constant. Similarly, \( v=y+f(x,y) \) cannot be constant, since otherwise \(G=q(v)\) would be constant and again \( \operatorname{eig} JH=0 \) would follow.

The outer polynomials must be affine

From \( {}^1p(u)=a \) and the non-constancy of \(u\), it follows that \({}^1p\) itself must be constant.

Indeed, the substitution homomorphism

is injective because \(u\) is non-constant. Hence

implies the polynomial identity

Therefore, for some constant \(\alpha\in\mathbb C\), the polynomial \(p\) has the form

Likewise, for some constant \(\beta\in\mathbb C\), \( {}^1q(v)=b \) implies

Since \(a,b\ne0\), both \(p\) and \(q\) are invertible affine polynomials.

In particular, the Keller condition excludes all genuinely non-linear outer polynomials in this class:

Differential equation for \(f\)

Equation (41) gives

To solve this polynomial differential equation, introduce the linear coordinates

and

Then

For the differential operators,

and

Therefore,

Equation (46) becomes

Integration with respect to \(s\) yields

where \( h\in\mathbb C[t] \) is an arbitrary polynomial.

In the original coordinates, the general solution is

Since

this can also be written as

Complete form of the Keller maps

Equations (43), (44), and (53) give

and

Thus every Keller map in the class under consideration has the form

and

where

and

Conversely, every map of the form (58)–(60) satisfies the Keller condition

Proof.

Equation (53) gives

and

Hence

Moreover, \( {}^1p=a \) and \( {}^1q=b. \) Equation (36) therefore gives

This also proves the converse.\(\square\)

Representation in sum and difference coordinates

Again set

For the inner coordinates

and

one first has

Furthermore, equation (53) gives

Thus, in the coordinates \((r,s)\), the inner map

has the particularly simple triangular form

The difference coordinate \(r=x-y\) remains unchanged, whereas the sum coordinate \(s=x+y\) is scaled by \(c\) and shifted by a polynomial in \(r\).

Since \(c\ne0\), this triangular map is polynomially invertible.

Explicit inverse map

Let

be the image coordinates.

Since

and

the outer affine maps are first inverted:

and

Then

and

Subtraction yields

Equation (62) gives

Hence

With

one obtains

Finally, equation (49) gives

and

Thus the explicit inverse map is

where \( u=\tilde a\,(U-\alpha) \) and \( v=\tilde b\,(V-\beta). \) Since \(h\) is polynomial and \(a,b,c\ne0\), the inverse \(H^{-1}\) is polynomial as well.

Direct exclusion of double points

Suppose that there exist two points

such that

Since \(p\) and \(q\) are affine polynomials with non-zero linear coefficients, equation (71) implies

and

Subtracting equations (72) and (73) gives

Thus the difference coordinates coincide: \( r_1=r_2. \) Consequently,

Adding equations (72) and (73) gives

Using equation (75), one obtains

Since \(c\ne0\),

Equality of the sum and difference coordinates now gives \( x_1=x_2 \) and \( y_1=y_2. \) This contradicts the assumption of a genuine double-point pair. Therefore \(H\) is injective.\(\square\)

The normalised case

In the normalised Keller case \( e=1, \) one has

Hence

The general form of \(f\) is

If, in addition, the outer polynomials are normalised by \( p(t)=t \) and \( q(t)=t, \) then \( a=b=1 \) and therefore \( c=1. \) In this case, \(f\) reduces to

The map is then

Its inverse map is

Indeed, \( F-G=x-y \) remains unchanged.

Shear representation

In the coordinates

the inner map

has the form

It is therefore a triangular map consisting of a scaling of the second coordinate and a shear:

The complete map \(H\) is then obtained by the affine transformations

in the two image components.

Thus every Keller map in this class is a composition of

and

Summary

Let \( f\in\mathbb C[x,y] \) and \( p,q\in\mathbb C[t]. \) Define

and

Then

holds if and only if there exist constants

and a polynomial \( h\in\mathbb C[t] \) such that

and

In this case, \(H\) has the explicit polynomial inverse (70). In particular, \(H\) is a polynomial automorphism and has no genuine double-point pair.

Within this comprehensive class, the Keller condition therefore forces both the affineness of the outer polynomials \(p\) and \(q\) and the complete sum–difference structure of the inner polynomial \(f\).\(\square\)

Remark: The case \(F(x,y)=x+f(x,y)\) and \(G(x,y)=y+p(F(x,y))\) is easier to prove.

Main theorem

Every two-component polynomial map \( H:\mathbb{C}^2\longrightarrow\mathbb{C}^2 \) of degree at most \(2\) satisfying \( \operatorname{eig} JH=e\ne0 \) is a polynomial automorphism. More precisely,

In particular, no genuine double-point pair exists up to degree \(2\).

In the quadratic case, the infinitesimal order calculation yields

Since a term of order \(\delta^2\) cannot cancel a non-zero term of order \(\delta\), it follows successively that \( \operatorname{tr}JR=0 \) and \( \operatorname{eig} JR=0. \) The resulting quadratic map is a shear and has the explicit polynomial inverse (25).\(\square\)

Scope of the proof

The preceding proof applies fully to \( \deg H\le2. \) For a map of higher degree, the order lemma separates the different homogeneous contributions to the identity \( \operatorname{eig} JH=e, \) but it does not by itself prove that all higher homogeneous parts of \(H\) vanish. An extension to arbitrary degrees therefore additionally requires either a valid reduction to the quadratic case or a separate exclusion of the higher homogeneous parts.\(\square\)