Linear Programming

The following likewise presupposes Set Theory, Topology and Nonstandard Analysis.

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\)

Theorem: For all regular \(C\;(= A^HA)\) with \(n\) as before, \(C^{-1} \in {}^{\nu}\mathbb{C}^{n\times n}\) can be computed via the Schur complement² and identity matrices \(I\) multiplied with \(\varepsilon \in \mathcal{O}(\tilde{\nu})\) in running time \(T(n) = \mathcal{O}(p^2)\), since \[ T(n) = \hat{T}(\check{n}) + \mathcal{O}(p^2) = 4T(\tilde{4}n) +\mathcal{O}(p^2 + \tilde{2}p^2) = \dots = \mathcal{O}(p^2) \] holds by virtue of the GS.\(\square\)

Corollary: Every solvable linear system \(Ax = b\) can be solved for regular or, as above, \(\varepsilon\)-perturbed singular \(A^HA\) with accuracy in \(]0, 1[\) in running time \(\mathcal{O}(p^2)\).\(\square\)

Theorem: The intex method solves every solvable LP with inter-/extrapolations in \(\mathcal{O}({}_2(\tilde{\alpha}\rho)^2 dmn)\).

Proof and algorithm: Let \(z := m + n\) and \(d \in [0, 1]\) be the density of \(A\). First we normalise and scale \[ b^{T}y – c^{T}x \le 0,\quad Ax \le b,\quad A^{T}y \ge c. \] Let \[ P_r := \{(x, y)^T \in {}^{\nu}\mathbb{R}_{\ge 0}^{z} : b^{T}y – c^{T}x \le r \in [0, \rho],\; Ax – b \le r{\upharpoonleft}_m,\; c – A^{T}y \le r{\upharpoonleft}_n\} \] have radius \[\rho := s|{\min} \{b_1, …, b_m, -c_1, …, -c_n\}|\] and scaling factor \(s \in [1, 2]\). Then \(0{\upharpoonleft}_z \in \partial P_{\rho}\) holds. By the strong duality theorem³ the LP \[ \min \{ r \in [0, \rho] : (x, y)^T \in P_r\} \] with the computational accuracy \(\alpha\) also solves the LPs \[ \max \{{c}^{T}x : x \in {}^{\nu}\mathbb{R}_{\ge 0}^{n}, Ax \le b\} \quad\text{and}\quad \min \{{b}^{T}y : y \in {}^{\nu}\mathbb{R}_{\ge 0}^{m}, A^{T}y \ge c\}. \] Its solution is the geometric centroid \(g\) of the polytope \(P_0\). With \[ v_k^* := \min \check{v}_k + \max \check{v}_k,\quad k = 1, …, \grave{z}, \] \(g\) is approximated by \(v_0 := (x_0, y_0, r_0)^T\) until \(||\Delta v||_1\) is sufficiently small. The solution \(t^o(x^o, y^o, r^o)^T\) of the two-dimensional LP \[ \min \{r \in [0, \rho] : t \in {}^{\nu}\mathbb{R}_{> 0},\; t(x_0, y_0)^T \in P_r\} \] approximates \(g\) better. After obtaining \(r \le \check{\rho}\), repeat the procedure with \(t^o(x^o, y^o)^T\) until, if applicable, \(g \in P_0\) has been computed in \(\mathcal{O}({}_2(\tilde{\alpha}\rho)^2 dmn)\). Solving all two-dimensional LPs \(\min_k r_k\) by bisection for \(r_k \in {}^{\nu}\mathbb{R}_{\ge 0}\), \(k = 1, …, z\), each in \(\mathcal{O}({}_2\tilde{\alpha}dmn)\), yields \(u \in {}^{\nu}\mathbb{R}^k\) with \( u_k := \Delta v_k \Delta r_k/r\) and \(r := \min_k \Delta r_k. \) For simplicity assume \(|\Delta v_1| = \dots = |\Delta v_{z}|\). Solve likewise in this setting \(\min r_{\grave{z}}\) for \(v^* := v + wu\) with \(w \in {}^{\nu}\mathbb{R}_{\ge 0}\). If \(\min_k \Delta r_k\, r = 0\) follows, stop computing; otherwise repeat until \(\min r = 0\) or \(\min r > 0\) is established.\(\square\)

Remark: Deliberate rounding errors in \(b\) and \(c\) can make a solution unique. If the method is optimised for distributed computing in \({}^{\nu}\mathbb{R}^{\nu}\), the running time is only \(\mathcal{O}(1)\). If \(\acute{q}_0 \le \lfloor q_j \rfloor – q_j\) holds for integer-valued solutions \(q_j\), the method performs well for (mixed-)integer problems and (non-)convex (Pareto) optimisation. If necessary, constraints can temporarily be relaxed. The following statements can also be transferred to the complex case without difficulty:

A primal-dual merit method for LPs

Let \(A\in{}^{\nu}\mathbb{R}^{m\times n}\), \(b\in{}^{\nu}\mathbb{R}^m\), \(c\in{}^{\nu}\mathbb{R}^n\), and \(e\in{}^{\nu}\mathbb{R}\). Consider the LP \[ \max \; e + c^\top x \qquad\text{subject to}\qquad Ax \le b,\qquad x \ge 0. \] A dual vector \(y\ge 0\), representing the dual condition \(A^\top y \ge c\), is introduced alongside the primal vector \(x\ge 0\). The method works directly in the primal-dual space \((x,y)\).

In contrast to classical simplex or interior-point methods, the objective function is not maximised directly. Instead, a merit function is first minimised, jointly capturing primal violations, dual violations, and the duality gap. The homogenised auxiliary component remains fixed at \(z=1\) and is not treated as a free variable. In practical implementation, two closely related variants are used.

The standard variant employs a balanced diagonal scaling, a scaled primal-dual merit function, and a two-phase procedure with symmetric final polish. It solves the majority of LPs considered in practice. For more difficult or numerically sensitive cases, an extended variant is used which retains the basic idea of the merit method but incorporates additional complementarity terms, adaptive weights, blockwise updates, and multi-stage refinement steps.

Primal-dual residuals

For \(x\ge 0\) and \(y\ge 0\), define the primal and dual residuals by \[ r^p(x):=\max(Ax-b,0)\in{}^{\nu}\mathbb{R}^m, \qquad r^d(y):=\max(c-A^\top y,0)\in{}^{\nu}\mathbb{R}^n, \] where the maximum is to be understood componentwise. In addition, the signed duality gap \[ g(x,y):=b^\top y-c^\top x \] is used; in fact, only its positive violation is penalised: \[ g_+(x,y):=\max(g(x,y),0). \] Thus the unscaled error quantities are \[ p_\infty=\|r^p(x)\|_\infty,\qquad d_\infty=\|r^d(y)\|_\infty,\qquad g_+=g_+(x,y). \] Moreover, \[ p^*=c^\top x,\qquad d^*=b^\top y \] are the primal and dual objective values respectively; the constant offset \(e\) is added at the end.

In both implementation variants, these quantities form the central numerical control basis. In the standard variant they are monitored both in the transformed problem and in the original one; in the extended variant additional complementarity measures appear, capturing the mutual compatibility of primal slacks and dual variables, and of primal variables and dual slacks.

Balanced scaling

To improve stability, the standard variant does not work directly with \((x,y)\) in the first phase, but rather with scaled variables \[ x=S_x\breve x,\qquad y=S_y\breve y, \] where \[ S_x=\operatorname{diag}(s_x)\in{}^{\nu}\mathbb{R}^{n\times n}, \qquad S_y=\operatorname{diag}(s_y)\in{}^{\nu}\mathbb{R}^{m\times m} \] are diagonal scaling matrices.

In addition, diagonal residual weights \[ D_p=\operatorname{diag}(d_p)\in{}^{\nu}\mathbb{R}^{m\times m}, \qquad D_d=\operatorname{diag}(d_d)\in{}^{\nu}\mathbb{R}^{n\times n} \] and a scalar gap weight \(d_g>0\) are constructed. These quantities are determined by several balancing passes. Starting from \[ d_p=1{\upharpoonleft}_m,\qquad d_d=1{\upharpoonleft}_n,\qquad s_y=1{\upharpoonleft}_m,\qquad s_x=1{\upharpoonleft}_n \] the rows \[ d_{p,i}^{-2} = b_i^2+{\LARGE{\textbf{+}}}_{j=1}^n (A_{ij}s_{x,j})^2, \qquad i=1,\dots,m, \] \[ d_{d,j}^{-2} = c_j^2+{\LARGE{\textbf{+}}}_{i=1}^m (A_{ij}s_{y,i})^2, \qquad j=1,\dots,n, \] are repeatedly scaled, where \[ (\gamma\tilde d_g)^2 = {\LARGE{\textbf{+}}}_{i=1}^m (b_i s_{y,i})^2 + {\LARGE{\textbf{+}}}_{j=1}^n (c_j s_{x,j})^2 \] holds and is updated with the parameter \(\gamma=\texttt{gap_weight}\). This is followed by the variable scalings \[ s_{y,i}^{-2} = (d_g b_i)^2+{\LARGE{\textbf{+}}}_{j=1}^n (d_{d,j}A_{ij})^2, \qquad i=1,\dots,m, \] \[ s_{x,j}^{-2} = (d_g c_j)^2+{\LARGE{\textbf{+}}}_{i=1}^m (d_{p,i}A_{ij})^2, \qquad j=1,\dots,n. \] If a zero-norm case occurs, the corresponding scaling is set to \(1\).

This balanced scaling is characteristic of the standard variant. It typically improves the balance between the primal and dual residuals and provides a robust starting point for the subsequent unscaled refinement. By contrast, the extended variant uses a simpler direct normalisation of the error terms and does not employ this explicit diagonal balancing procedure.

The scaled merit function

In the scaled variables \(\breve x,\breve y\), the residuals are \[ r^p(\breve x)=\max(AS_x\breve x-b,0), \qquad r^d(\breve y)=\max(c-A^\top S_y\breve y,0), \] as well as \[ g(\breve x,\breve y)=b^\top S_y\breve y-c^\top S_x\breve x, \qquad g_+(\breve x,\breve y)=\max(g(\breve x,\breve y),0). \] The minimised merit function is \[ \hat{\Phi}(\breve x,\breve y) = \|D_p\,r^p(\breve x)\|_2^2 + \|D_d\,r^d(\breve y)\|_2^2 + d_g^2\,g_+(\breve x,\breve y)^2. \] In a slightly regularised implementation variant, an additional very small term \( \check{\varepsilon}\bigl(\|\breve x\|_2^2+\|\breve y\|_2^2\bigr) \) may be included without altering the basic structure of the method.

Thus, in its essential form, the merit function contains exactly three contributions:

• primal violation,
• dual violation,
• positive gap violation.

Complementarity terms are no longer penalised separately in the final version.

This representation describes the basic form of the standard variant. The extended variant uses the same primal-dual idea, but augments the merit function by additional quadratic complementarity terms. In particular, products of positive primal slacks with dual variables, and of primal variables with positive dual slacks, are penalised separately there. Moreover, in this variant the gap is treated quadratically in a symmetric manner, so that both negative and positive gap deviations enter the objective function equally. Nevertheless, the above formula remains the conceptual starting form of the simpler method, which is decisive for most problems.

Gradient of the merit function

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 with respect to the scaled variables 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. \] In the regularised variant, the additional summands \(\varepsilon\breve x_j\) and \(\varepsilon\breve y_i\) are included.

Hence the primal side is determined by the violations of \(Ax\le b\), and the dual side by the violations of \(A^\top y\ge c\); the gap term couples both blocks.

This gradient also describes the standard variant. In the extended variant, derivatives of the complementarity terms are added. As a result, that method is governed not only by primal, dual, and gap errors, but, once feasibility has improved sufficiently, also by the requirement that positive primal slacks should not occur simultaneously with positive dual variables, and positive primal variables should not occur simultaneously with positive dual slacks.

Projected gradient descent

Starting from \((\breve x^k,\breve y^k)\), a joint projected descent step is performed: \[ x^\prime = \Pi_{\mathbb{R}_+^n}\bigl(\breve x^k-\alpha_k\nabla_{\breve x}\hat{\Phi}(\breve x^k,\breve y^k)\bigr), \qquad y^\prime = \Pi_{\mathbb{R}_+^m}\bigl(\breve y^k-\alpha_k\nabla_{\breve y}\hat{\Phi}(\breve x^k,\breve y^k)\bigr), \] where \(\Pi_{\mathbb{R}_+^\ell}\) denotes the orthogonal projection onto the non-negative orthant.

The step is accepted only if an Armijo condition is satisfied: \[ \hat{\Phi}(x^\prime,y^\prime) \le \hat{\Phi}(\breve x^k,\breve y^k) – \tilde{\alpha}_k c_1\left( \|x^\prime -\breve x^k\|_2^2+\|y^\prime -\breve y^k\|_2^2\right), \] where \(c_1>0\) is the Armijo parameter. Otherwise, the step size is halved until the condition is met or a lower bound is reached.

This description applies in particular to the standard variant. In the extended variant, the same basic idea is realised blockwise: there, the \(x\)- and \(y\)-steps are treated separately, each projected, safeguarded by backtracking, and continued with their own step sizes. In addition, an optional momentum extrapolation may be used there, producing a Nesterov-like acceleration as long as the merit function does not increase again.

Barzilai-Borwein updating

After an accepted step, the step size is updated in the spirit of Barzilai-Borwein. For this purpose, \[ s_k= \begin{pmatrix} \breve x^{k+1}-\breve x^k\\[0.3em] \breve y^{k+1}-\breve y^k \end{pmatrix}, \qquad q_k= \begin{pmatrix} \nabla_{\breve x}\hat{\Phi}(\breve x^{k+1},\breve y^{k+1})-\nabla_{\breve x}\hat{\Phi}(\breve x^k,\breve y^k)\\[0.3em] \nabla_{\breve y}\hat{\Phi}(\breve x^{k+1},\breve y^{k+1})-\nabla_{\breve y}\hat{\Phi}(\breve x^k,\breve y^k) \end{pmatrix} \] are formed. If \(\eta_k := s_k^\top q_k>0\), then \( \alpha_{k+1} = \tilde{\eta}_k{s_k^\top s_k} \) is used and projected onto a prescribed interval \([\alpha_{\min},\alpha_{\max}]\). If \(\eta_k\le 0\), the last usable step size is retained.

In the standard variant, this updating is performed for the joint primal-dual step. In the extended variant, the same idea is applied blockwise to the individual \(x\)- and \(y\)-updates. In substance, however, the objective remains the same in both cases: to obtain a locally adaptive step size that is largely robust under scaling, without requiring an explicit Hessian matrix.

Multi-phase strategy, back-transformation, and termination

The method operates in three interrelated phases.

Phase 1.
First, the problem is solved with balanced scaling. The number of balancing passes is specified by the parameter balance_passes. This phase usually reduces the transformed residuals very reliably and provides a robust primal-dual approximation point.

Phase 2.
If Phase 1 has not yet produced a sufficiently accurate solution of the original problem and iterations remain, a second run is started from the physical warm start \( x=S_x\breve x\) and \(y=S_y\breve y \). This run is performed on the original problem without balanced scaling; it serves as a warm-start polish of the unscaled primal and dual residuals and of the duality gap. This avoids noticeable errors in the original model when the transformed model has already been solved well.

Phase 3.
If Phase 2 still does not provide an exact LP solution, a basis polish follows: the computed approximation point generates candidates for active constraints \(a_i^\top x=b_i\) and non-negativity conditions \(x_j=0\). From these, linearly independent basic configurations of size \(n\) are formed and tested, and an LS for \(x\) is solved. A consistent dual vector \(y\) is then reconstructed. The best basic solution is retained, which functionally corresponds to a crossover from the smooth merit phase to an exact basic solution.

This threefold division describes the standard variant to a suitable approximation. There, the third step consists in practice of a symmetric final polish which extracts a more consistent active substructure from the pattern of small primal and dual slacks and positive variables, and reconstructs an improved final solution from it. The term basis polish is therefore to be understood functionally, not in a simplex-theoretical sense. By contrast, the extended variant does not pursue any explicit basis reconstruction.

After an initial feasibility and gap phase, several adaptive refinement stages are carried out there, in which the weights of the primal, dual, and gap terms, as well as the strength of the complementarity penalties, are adjusted step by step. Short repair runs without complementarity penalty may stabilise feasibility if an overly strong emphasis on complementarity temporarily worsens primal or dual error.

Back-transformation.
After each scaled phase, the physical variables are recovered by \( x=S_x\breve x\) and \( y=S_y\breve y \). In the unscaled second phase this back-transformation is trivial.

Termination.
In Phase 1 and Phase 2, the method terminates if either \( \max\{p_\infty,d_\infty,g_+\}\le \tau \) holds, or the merit function has become sufficiently small, namely \( \hat{\Phi}(\breve x,\breve y)\le \tau^2. \) If no exact feasible solution has been reached by then, Phase 3 determines the final basic solution. In the extended variant, termination is analogous via small primal, dual, and gap errors; in addition, complementarity measures are monitored there.

Merit is thus overall a multi-phase first-order primal–dual procedure with projection, backtracking, Barzilai–Borwein step-size control, and a algebraic end polish or adaptive refinement. Its strength lies in the simple structure of the first two phases, in the direct control of primal and dual violation as well as of the duality gap, and in the possibility of reconstructing a highly accurate final solution from a smooth approximation point.

Complexity and assessment

The method solves LPs in a non-standard way – unlike simplex or interior-point methods. Methodologically, it is simple, transparent, and readily modifiable. Each evaluation of the merit function and its gradient is determined by a constant number of matrix-vector-like passes over \(A\) and \(A^\top\). All three phases are essential for practical robustness. For sparse matrices, the cost is proportional to the number of non-zeros, i.e.\ \(\mathcal{O}(\mathrm{nnz}(A))\), and for dense representation it is \(\mathcal{O}(mn)\).

Let \(K\) be the number of accepted iterations, \(B\) the total number of backtracking evaluations, and \(P\) the number of balancing passes. Then the total runtime of the first two phases is of the order of a multiple of \( (K+B+P)\,\mathcal{O}(\mathrm{nnz}(A)) \) or, in dense representation, \( (K+B+P)\,\mathcal{O}(mn). \) The additional memory requirement beyond storing \(A\) is only linear in the dimensions in these phases, namely \(\mathcal{O}(m+n)\); together with matrix storage this gives \(\mathcal{O}(mn+m+n)=\mathcal{O}(mn)\) in the dense case.

The method works independently of pivot rules. At the end, a basis polish generates from a small set of suspected active constraints and non-negativity conditions a finite number of basis candidates. For each candidate, an LS of size \(n\) has to be solved, which costs \(\mathcal{O}(n^3)\) in dense treatment. If \(Q\) candidate equations are considered and \(T\) basis combinations are actually tested, then the third phase contributes approximately an additional cost of \(\mathcal{O}(T\,n^3)\), with \(T\le \binom{Q}{n}.\)

For the standard variant, the linear algebra of the first two merit phases dominates initially, before the symmetric final polish on a small active substructure causes additional, but typically much smaller, costs. The small, heuristically selected candidate set limits the number of tested combinations and is suitable for small and medium values of \(n\).

For the extended variant, by contrast, the explicit basis reconstruction step is omitted; there the additional computation time is determined by several adaptive stages, separate \(x\)- and \(y\)-updates, complementarity tests, and, where necessary, repair runs. Each iteration remains of matrix-vector-like complexity and makes the method simpler than those that solve large-dimensional LSs at every step.

Overall, the method replaces basis changes and Newton steps by the minimisation of a primal-dual residual function, and subsequently supplements this either by a targeted algebraic end polish or by adaptive complementarity-oriented refinements. Its efficiency depends more strongly on step sizes, scaling, and parameter choice than that of the established standard methods. The extended variant serves primarily for more difficult or numerically sensitive instances.

If an additional algebraic structure is present (Toeplitz, Hankel, block-circulant, or more generally convolution-like), then a geometric search is possible, but typically less efficient than a direct structured factorisation. In that case, the FTD and FSR methods prove complementary to the merit method, which is particularly suitable for problems with weak or absent algebraic structure.

Corollary: The LP     \( \max \{{||x – {x}^{o}||}_1 : {c}^{T}x = {c}^{T}{x}^{o},\; Ax \le b,\; x – {x}^{o} \in {[-1, 1]}^{n},\; x \in {}^{\nu}\mathbb{R}_{\ge 0}^{n}\} \)     can, given the first solution \({x}^{o}\), determine a second one in \(\mathcal{O}(mn)\), if applicable, where \({y}^{o}\) can be treated analogously.\(\square\)

Corollary: The LP    \( \max \{\nu|\lambda_j| + ||x_j||_1 : Ax_j = \lambda_j x_j,\; |\lambda_j| \in [0, r_j],\; x_j \in {[-1, 1]}^{n*}\} \)\qquad can, for each \(j = 1, …, n\), determine the eigenvalue \(\lambda_j \in {}^{\nu}\mathbb{R}\) for the eigenvector \(x_j \in {}^{\nu}\mathbb{R}^{n}\) of \(A \in {}^{\nu}\mathbb{R}^{n \times n}\) in \(\mathcal{O}(n^2)\).\(\square\)

Corollary: The LPs \(\max \{x_{j} : Ax = 0\}\) yield all solutions of \(Ax = b\) in \(\mathcal{O}(mn)\). The matrix \(A\) is regular if and only if the LP \(\max \{||x||_{1} : Ax = 0\} = 0\).\(\square\)

Corollary: Let \(\alpha_{j} := A_{.j}^{-1}\) for \(j = 1, …, n\) with respect to the matrix \(A^{-1} \in {}^{\nu}\mathbb{R}^{n \times n}\), and let \(\delta_{ij}\) denote the Kronecker delta. For regular \(A\) and eigenproduct \(\neq 0\), each \( A \alpha_{j} = (\delta_{1j}, …, \delta_{nj})^{T} \) can be solved in \(\mathcal{O}(n^2)\).\(\square\)

Corollary: The polynomial \(c_{\acute{n}}x^{\acute{n}} + \dots + c_1x + c_0\) can be determined in \(\mathcal{O}(n^2)\) from measurements \((x_j, y_j) \in {}^{\nu}\mathbb{R}^2\) and residuals \(r_j \in {}^{\nu}\mathbb{R}\) for \(j = 1, …, m\), as well as \((c, x)^T \in {}^{\nu}\mathbb{R}^{\grave{n}}\) (also under additional side constraints), from \[ y = r + Xc,\quad X^Tr = 0,\quad X = (x_j^k)_{j,k},\ k = 0, …, \acute{n}.\square \] Corollary: The DFT (see Scientific Computing) can determine the Padé approximation

\(\frac{a(x)}{b(x)} = {\LARGE{\textbf{+}}}_{m=0}^n a_m x^m / {\LARGE{\textbf{+}}}_{m=0}^n b_m x^m\) of \(f(x)\) for \(x \in {}^{\nu}\mathbb{R}\) in \(\mathcal{O}(n^2)\).\(\square\)

Corollary: Merit methods and two-dimensional bisection or Newton methods solve most solvable convex programmes \[ \min \{f_{1}(x) : x \in {}^{\nu}\mathbb{R}^{n},\ (f_{2}(x), …, f_{m}(x))^{T} \le 0\} \] with convex functions \(f_{k} \in {}^{\nu}\mathbb{R}\) for \(k = 1, …, m\) in polynomial time if the number of operands \(x_{j}\) of the \(f_{k}\) is \(\le {\nu}^{\nu-3}\) and there exists an \(x\) with \(f_{k}(x) < 0\) for all \(k > 1.\square\)

Corollary: The LP     \( \max \{x : y_{\acute{m}} – xy_m = b_{\acute{m}},\; y_n = y_0 = 0,\; x \le x_0 \in {}^{\nu}\mathbb{R}\} \)     can, for \(m = 1, …, n\), determine an \(b_{\acute{n}}x^{\acute{n}} + \dots + b_1x + b_0 = 0\) solving \(x \in {}^{\nu}\mathbb{R}\) according to Horner’s scheme in \(\mathcal{O}(n)\).\(\square\)

The FTD is based on the observation that a regular matrix \(A\in{}^{\nu}\mathbb{R}^{n\times n}\), whose entries vary only slowly along the diagonals or can be described by convolution operators, can be well approximated as the product of two Toeplitz matrices \(T_{1},T_{2}\). At the beginning, \(x_{0}=\max|A_{ij}|\) is chosen; subsequently, the sequences \(x_{k},y_{k}\) are formed by local arithmetic means of neighbouring entries. The diagonal sums \[ z_{1}=x_{0}y_{0}+x_{1}y_{1}+\ldots,\qquad z_{2}=x_{0}y_{-1}+x_{1}y_{0}+\ldots \] are computed recursively until \(A\approx T_{1}T_{2}\) is reached. The normalisation of the columns of \(T_{1}\) and \(T_{2}\) stabilises the recursion. By backward substitution along the diagonals, a unique approximation is obtained. The number of steps per diagonal is \(\mathcal{O}(\ell)\). If \(T_{1}^{-1}\) and \(T_{2}^{-1}\) are then computed, the purely structured operations imply \[ T(n)=\mathcal{O}(p^{2}),\qquad p=\ell n. \] For symmetric systems, the factorisation \(H\approx P^{T}DP\) with diagonal matrix \(D\) and block-Toeplitz matrix \(P\) arises from \(H=T_{1}^{T}T_{1}\), and is used in the LUX method for solving structured LPs. It replaces classical Cholesky or QR methods and preserves Toeplitz, Hankel, and circulant structures exactly.

In parallel, the FSR methods are based on recursive shifts and modulo reconstruction. A matrix \(A\) is decomposed into low- and high-weight parts, \[ A=(A_{\mathrm{hi}}\ll s)+A_{\mathrm{lo}}, \] where \(\ll s\) denotes a generative shift. The products of the individual parts are computed recursively and then recombined modulo a suitably chosen \(s\) into an exact result. In this way, all intermediate values remain dynamically small, the recursion is fully parallelisable, and BigInt and BigFloat arithmetic do not lead to overflow problems. Thus, for matrix multiplication, inversion, and Schur complement steps, the same complexity holds as for FTD, \[ T_{\mathrm{FSR}}(n)=\mathcal{O}(p^{2}). \] The three methods complement one another naturally. FTD solves structured LSs and symmetric factorisations exactly and requires only convolution operators and diagonal updates. FSR performs the arithmetic operations recursively and in a memory-saving fashion and is suitable for very large integers and high precision. If a Toeplitz, Hankel, or general convolution structure is present in an LP or an LS, then the composition \( A\approx T_{1}T_{2},\qquad T_{1}^{-1},T_{2}^{-1}\) recursively by FSR, projection or LP section by merit, reduces the overall effort to \(\mathcal{O}(p^{2})\) for the structured LSs and the LPs. The methods work independently of the word length of the entries and scale almost linearly for numbers of size \(\mathcal{O}(\nu)\).

The combination of FTD, FSR, and merit thus forms a universal toolbox for structured linear algebra and optimisation. It combines structure (Toeplitz, Hankel, circulant), recursion (shift, convolution), and geometry (projection) into a procedure that guarantees both regularity and exact solutions in short runtime. It therefore offers a unified approach to problems that have hitherto been solvable only separately by numerical or geometric methods.

The following are orders of magnitude of the new methods with \(p=\ell n\) and newly computable ranges. The stated thresholds are to be understood heuristically and depend on implementation, hardware, and the desired accuracy; however, they indicate the ranges in which the new methods could systematically outperform the classical approaches.

Please, input real linear programme (separators spaces, last row objective function, first column right-hand sides as in the example):

code of the simplex method

© 2008-2026 by Boris Haase

top

1. see Golub, Gene H.; van Loan, Charles F.: Matrix Computations; 3rd Ed.; 1996; Johns Hopkins University Press; Baltimore, p. 31 ff.
2. see Knabner, Peter; Barth, Wolf: Lineare Algebra; 2., überarb. und erw. Aufl.; 2018; Springer; Berlin, p. 223
3. Vanderbei, Robert J.: Linear Programming; 3rd Ed.; 2008; Springer; New York., p. 60 – 65