Solving Non-Symmetric Linear Systems: A Guide to GMRES and BiCGSTAB

When faced with linear equations of the form $A\mathbf{x} = \mathbf{b}$, mathematicians and engineers often rely on a wealth of powerful tools. However, many of the most efficient and elegant methods are designed for matrices with a special property: symmetry. The moment this symmetry is lost, we enter a more complex and challenging realm. This article addresses the critical question: how do we effectively solve large linear systems when the underlying matrix $A$ is non-symmetric? The breakdown of standard techniques necessitates new strategies, each with a unique philosophy for navigating the problem's landscape.

This guide provides a conceptual journey into the world of iterative solvers for non-symmetric systems. We will move beyond simply trying to force symmetry and instead explore methods designed to work directly with the problem's inherent structure. Across the following chapters, you will discover the core mechanics of two celebrated algorithms and see how they connect to real-world phenomena. In "Principles and Mechanisms", we will delve into the inner workings of the robust GMRES and the nimble BiCGSTAB methods, contrasting their perfectionist and pragmatist approaches to finding a solution within Krylov subspaces. Following that, "Applications and Interdisciplinary Connections" will reveal where these non-symmetric problems arise, from the flow of heat and pollutants in engineering to the fundamental mechanics of materials and the structure of economic networks.

Now that we have been introduced to the vast stage of non-symmetric linear systems, let’s peel back the curtain and look at the gears and levers that make the machinery work. How do we actually begin to solve an equation like $A\mathbf{x} = \mathbf{b}$ when the beautiful, orderly world of symmetric matrices is no longer our home? The story of these methods is a fascinating journey of brilliant ideas, frustrating pitfalls, and clever compromises.

The first and most natural impulse, when faced with an unfamiliar and difficult problem, is to try to turn it into a familiar and easy one. We have powerful, efficient tools like the Conjugate Gradient (CG) method for matrices that are ​​symmetric and positive-definite​​. So, can't we just force our non-symmetric matrix $A$ to become symmetric?

Indeed, we can! Consider our original equation, $A\mathbf{x} = \mathbf{b}$. If we simply multiply both sides by the transpose of our matrix, $A^T$, we get a new equation:

Let's look at the new matrix in charge, $A^T A$. It is a mathematical fact that for any invertible matrix $A$, the product $A^T A$ is always symmetric and positive-definite. Suddenly, we are back in familiar territory! We can unleash the powerful CG method on this new system to find our solution $\mathbf{x}$. This "trick," known as the method of ​​normal equations​​, seems like a perfect, elegant solution. It's one of a couple of ways you can transform the problem to regain symmetry.

But here lies a subtle and dangerous trap. While we've made the problem look easier, we've often made the underlying numerical challenge much, much harder. The "difficulty" of a linear system is often measured by something called its ​​condition number​​. A high condition number means the system is sensitive and ill-behaved, like trying to weigh a feather on a windy day. The condition number of our new matrix, $A^T A$, is the square of the condition number of the original matrix $A$. If $A$ was already a bit tricky, $A^T A$ can be monstrously so. By forcing symmetry, we've amplified the very numerical difficulties we were trying to overcome. We need a more direct, more native approach.

Instead of contorting the problem, let's explore it. Imagine you have an initial guess for the solution, $\mathbf{x}_0$. It's probably wrong. The error, or ​​residual​​, is $\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0$. This residual vector points in a direction we need to go to correct our guess.

Now, what if we see where the matrix $A$ takes this direction? We get a new vector, $A\mathbf{r}_0$. And what if we do it again? We get $A^2\mathbf{r}_0$. We can think of this as a series of exploration steps. The initial residual $\mathbf{r}_0$ is our first step. Applying $A$ is like taking a second, transformed step based on the first, and so on.

The set of all locations we can reach by combining these steps, $\text{span}\{\mathbf{r}_0, A\mathbf{r}_0, A^2\mathbf{r}_0, \dots, A^{k-1}\mathbf{r}_0\}$, forms a mathematical playground called the ​​Krylov subspace​​. This subspace is the heart of modern iterative methods. The grand strategy is no longer to transform the matrix, but to find the best possible approximation to our true solution that lies within this expanding subspace. The various methods we'll discuss are simply different philosophies about what "best" truly means.

Our first philosophy is that of an uncompromising perfectionist: the ​​Generalized Minimal Residual Method (GMRES)​​. At every single iteration $k$, GMRES asks a simple, powerful question: "Of all the possible solutions in the $k$-dimensional Krylov subspace, which one makes the remaining residual $\mathbf{r}_k = \mathbf{b} - A\mathbf{x}_k$ as small as absolutely possible?" It seeks to minimize the length (the 2-norm) of this residual vector.

This philosophy has a beautiful and direct consequence. Since the Krylov subspace at step $k+1$ contains the entire subspace from step $k$, the search space for the best solution is always growing. Minimizing the same quantity over a larger space means the minimum can only get smaller or stay the same. Therefore, the residual norm in GMRES is ​​guaranteed to be monotonically non-increasing​​. The error will never get worse, which is wonderfully reassuring.

But how does it achieve this perfection? The engine inside GMRES is the ​​Arnoldi iteration​​. Think of it as an immaculate procedure that takes the raw, messy vectors $\{ \mathbf{r}_0, A\mathbf{r}_0, \dots \}$ that define our Krylov subspace and builds a pristine, orthonormal basis from them—a set of perfectly perpendicular unit vectors that span the same space. In doing so, it simultaneously builds a small, well-structured matrix (an upper Hessenberg matrix, $\tilde{H}_k$). The entire, high-dimensional problem of minimizing the residual is transformed into an equivalent, tiny least-squares problem involving this small matrix.

This perfectionism, however, comes at a cost.

• ​​Memory and Work​​: To find the optimal solution at step $k$, Arnoldi's process requires GMRES to keep all $k$ of the basis vectors it has so carefully constructed. For large systems where thousands of iterations might be needed, this can exhaust a computer's memory. This leads to a pragmatic compromise: ​​restarted GMRES​​, or ​​GMRES(m)​​, where the process is run for $m$ steps, and then completely restarted using the latest solution as a new guess.
• ​​The Perils of Restarts​​: This restart can be treacherous. By throwing away the basis vectors, we might be throwing away vital information about the landscape of the problem. For some tricky non-normal matrices, restarted GMRES can stagnate or fail to converge entirely. Imagine trying to navigate a complex maze but only being allowed to remember your last two turns. You could easily get stuck in a simple loop without ever reaching the exit. This exact scenario can be constructed, where GMRES(2) cycles forever, making no progress, while a full GMRES would have found the solution.
• ​​Stagnation​​: Even without restarts, GMRES can sometimes exhibit long periods of "stagnation." This often happens with non-normal matrices. The residual norm might decrease by an infinitesimal amount for dozens or hundreds of iterations before suddenly plummeting towards the solution. It is as if the algorithm is trudging across a vast, nearly flat plateau before finally finding the steep cliff that leads to the answer. For certain matrices, the reduction in error after the first step can be almost non-existent.

If GMRES is the perfectionist, our second philosophy is that of the resourceful pragmatist: the ​​Biconjugate Gradient Stabilized (BiCGSTAB)​​ method. This method abandons the expensive quest for a perfectly minimal residual at every step. In exchange, it gains tremendous efficiency in both memory and computation.

The story starts with its predecessor, BiCG. It ingeniously generalizes the Conjugate Gradient method by working with two sets of vectors instead of one. It introduces a "shadow" system involving the matrix transpose, $A^T$, and a ​​shadow residual​​, $\hat{\mathbf{r}}_0$. Instead of forcing a sequence of residuals to be orthogonal to each other (which only works for symmetric matrices), it forces the residuals $\mathbf{r}_k$ and the shadow residuals $\mathbf{r}^*_k$ to be "bi-orthogonal." The choice of the initial shadow residual is crucial; it sets the stage for the entire calculation and directly influences the path the algorithm takes. The purpose of this clever, two-sided bookkeeping is that it allows the algorithm to be built with "short recurrences"—each new direction depends only on the previous one, not the entire history. This is why it's so light on memory.

This cleverness, however, makes BiCG fragile. The bi-orthogonality condition can suddenly fail, leading to a division by zero—a catastrophic ​​breakdown​​. Even when it doesn't break down, its convergence can be wildly erratic, with the residual norm spiking up and down unpredictably.

This is where the "STAB" (stabilized) part saves the day. BiCGSTAB is a brilliant hybrid. As revealed in its very structure, each iteration is a two-act play:

1. ​​A BiCG Step​​: It first takes a provisional step using the efficient but potentially wild BiCG logic. Its direction is determined by the bi-orthogonality principle, governed by the shadow residual.
2. ​​A GMRES(1) Step​​: It then immediately "stabilizes" this move. It takes the result from the first step and performs a tiny, one-dimensional search to find the best correction along a new direction. This is, in essence, the smallest possible GMRES step imaginable.

This stabilization step acts as a damper, smoothing out the violent oscillations of pure BiCG. The convergence is no longer guaranteed to be monotonic—the residual norm can, and often does, increase temporarily—but the overall behavior is typically much faster and smoother on its way to the solution.

So we are left with two profoundly different, yet powerful, strategies for navigating the world of non-symmetric systems.

On one hand, we have ​​GMRES​​, the perfectionist. It is robust, its progress is guaranteed, but it demands a heavy price in memory and computation. Its restarted variant, GMRES(m), is a practical compromise, trading some of its robustness for feasibility, though one must be wary of stagnation and the pitfalls of a limited view.

On the other hand, we have ​​BiCGSTAB​​, the pragmatist. It is nimble, memory-efficient, and often very fast. It achieves this by giving up the guarantee of monotonic convergence, accepting a more unpredictable path to the solution in exchange for speed. It embraces a hybrid nature, combining the efficiency of one idea with the stability of another.

Which is better? There is no single answer. The choice is an art, guided by the specific properties of the matrix $A$, the computational resources at hand, and the desired accuracy. This journey into principles and mechanisms reveals the true beauty of numerical science: it is not just about computing an answer, but about the invention of elegant and diverse strategies, each with its own philosophy, strengths, and weaknesses, for finding a path through complexity.

In our previous discussion, we opened the physicist’s toolbox and examined a special class of tools: iterative solvers for non-symmetric linear systems. We familiarized ourselves with the elegant machinery of methods like GMRES and BiCGSTAB, designed to tackle equations where the underlying operator lacks the tidy, mirror-like quality of symmetry. It is a fair question to ask: where in the vast landscape of science and engineering do we encounter such mathematical beasts? Is this a niche pathology, a curious corner of mathematics, or does non-symmetry lie at the heart of phenomena we observe every day?

The answer, you might be delighted to discover, is that the ghost of non-symmetry is everywhere. While physics often celebrates symmetry—in the beautiful balance of conserved quantities, in the structure of crystals, in the action principles that govern motion—the moment we look at the real, bustling, and often messy world, we find direction. We find irreversibility. We find flow. A river flows to the sea; heat flows from hot to cold; influence spreads through a social network. These are one-way streets. And it is in the modeling of these directed processes that non-symmetric systems emerge not as an exception, but as the rule. Let us now go on a tour and see where they appear.

Our first stop is the most intuitive of all: the world of fluids, flows, and transport. Imagine releasing a drop of dye into a river. Two things happen. The dye spreads out in a growing, fuzzy cloud—this is ​​diffusion​​. It’s a beautifully symmetric process; the dye spreads equally in all directions from its center, driven by random molecular collisions. If the river were perfectly still, the problem of modeling this spread would give us a lovely symmetric matrix.

But the river is not still; it is flowing. The entire cloud of dye is also carried downstream. This is ​​convection​​ (or ​​advection​​). It is an inherently directed process. It gives the dye a collective “push” in a single direction. When we write down the equations of motion for our dye and translate them into a linear system for a computer to solve, the symmetric diffusion part is joined by a new, non-symmetric advection part. The total system is now non-symmetric. The matrix that describes the dye’s fate is no longer a mirror image of itself.

This simple picture of dye in a river scales up to grand challenges in science and engineering. The same advection-diffusion-reaction equations model the transport of pollutants in the atmosphere, the dispersal of nutrients in the ocean, and the flow of heat in industrial processes. In all these cases, the presence of a bulk flow—a wind, an ocean current, a forced coolant flow—introduces non-symmetry.

What’s more, the degree of non-symmetry has profound physical and numerical consequences. We can define a dimensionless quantity, the Péclet number, $Pe$, which measures the strength of convection relative to diffusion. When $Pe$ is small, diffusion dominates. The system is “mostly symmetric,” and our specialized solvers might not even be necessary. But when $Pe$ is large—in a strong wind or a fast-flowing river—convection dominates. The system becomes aggressively non-symmetric. Something more subtle happens, too: the system matrix can become highly non-normal. This is a particular kind of mathematical nastiness where the eigenvectors are far from orthogonal, a situation that is notoriously difficult for some of our elegant solvers. For instance, the otherwise robust GMRES method can see its convergence stagnate for many iterations, as if lost in a mathematical swamp, especially when it is “restarted” to save memory. In these very situations, the more erratic but persistent BiCGSTAB method can sometimes find its way to a solution more effectively.

The lesson here is profound: the physics dictates the mathematics. To solve a problem dominated by flow, we must use tools that respect the directed, non-symmetric nature of that flow. Even our choice of preconditioner—a helper matrix used to speed up convergence—must be physically informed. If we use a preconditioner that only understands the symmetric physics of diffusion, it will be utterly ineffective when faced with a problem dominated by non-symmetric convection.

Let us leave the world of fluids and turn to the seemingly static and stable world of solids. Surely a block of steel or the frame of a bridge can be described by symmetry? Often, yes. The analysis of complex structures, however, especially when they undergo large deformations, is a nonlinear problem. We solve it step-by-step using methods like the Newton-Raphson scheme, where each step requires solving a linear system. The matrix of this system, called the ​​tangent stiffness matrix​​ $\mathbf{K}_T$, tells us how the structure resists a small change in shape. And its symmetry, it turns out, is a direct reflection of the physics of the forces and the material itself.

First, consider the forces. A “dead load,” like gravity, is conservative. Its potential energy depends only on position, not the path taken. This conservatism leads to a symmetric tangent stiffness matrix. But what about a force whose direction depends on the object's orientation? Imagine the thrust from a rocket engine mounted on a flexible wing. As the wing bends, the direction of the thrust changes with it. This is a ​​follower force​​. Because it is not derivable from a simple scalar potential—the work it does depends on the path of the deformation—it breaks the variational structure of the problem. This physical non-conservatism manifests directly as a non-symmetric tangent stiffness matrix $\mathbf{K}_T$. To simulate such a system, we must wheel out our non-symmetric solvers.

Second, non-symmetry can arise from the deep, internal constitution of matter itself. When a metal is stressed beyond its elastic limit, it deforms permanently—this is ​​plasticity​​. In many models for materials like soils, rocks, and concrete, the rule that governs the onset of plastic flow (the yield function, $f$) is different from the rule that governs the direction of that flow (the plastic potential, $g$). This is called ​​non-associative plasticity​​. This subtle mismatch, deep within the constitutive law of the material, breaks the internal symmetry of its response. When we derive the consistent tangent matrix for the material, it emerges as non-symmetric. Thus, accurately predicting the behavior of a building during an earthquake or the stability of a tunnel wall requires us to solve non-symmetric linear systems, not because of the external forces, but because the very earth it is built on is non-symmetric in its response to stress!

The reach of non-symmetric systems extends far beyond continuum mechanics. Consider the discrete, interconnected world of networks. Think of the World Wide Web, an economic supply chain, or a food web in an ecosystem. These are networks of directed links: page A links to page B; company X supplies parts to company Y; the fox eats the rabbit.

Suppose we want to find a steady-state distribution of some quantity on such a network—the importance of a webpage, the economic output of a sector, the biomass of a species. This equilibrium is a balance between what flows in and what flows out. The analysis often leads to a classic linear system of the form $(\mathbf{I} - \mathbf{P})\mathbf{x} = \mathbf{s}$, where $\mathbf{P}$ is a transfer matrix describing how "stuff" is passed from node to node along directed edges, and $\mathbf{s}$ is an external source term. Because the underlying graph is directed, the matrix $\mathbf{P}$, and thus the system matrix $(\mathbf{I} - \mathbf{P})$, is inherently non-symmetric. The same mathematical structure can describe the equilibrium of an economy (where it is known as a Leontief input-output model) and the flow of nutrients in an ecosystem, revealing a surprising unity in the mathematics of interconnectedness.

Now for a final, crucial lesson in computational wisdom. The most famous non-symmetric system of our time is arguably the one that powers Google's PageRank algorithm. It seeks to find the "importance" vector $\mathbf{x}$ of all pages on the web, and it can be written precisely in the form $(\mathbf{I} - \alpha \mathbf{P})\mathbf{x} = \mathbf{b}$. Given that this is a non-symmetric system of staggering size—billions of equations—shouldn't we use our most sophisticated tool, like BiCGSTAB?

The answer, surprisingly, is no. Practitioners use a much simpler algorithm called the ​​power method​​. Why? Because while the problem is a non-symmetric linear system, it has a very special character. The power method, though slow to converge, is incredibly simple, robust, and requires minimal memory. Most importantly, it naturally preserves the physical meaning of the solution: every iterate remains a probability vector (all entries are non-negative and sum to one). BiCGSTAB, by contrast, would be plagued by the system's ill-conditioning, would require more memory and communication, and its intermediate solutions would be physically meaningless combinations with negative probabilities. The PageRank story is a beautiful reminder that true mastery is not just about knowing how to use powerful tools, but about understanding the specific character of a problem and knowing when a simpler tool is the right one for the job.

Our journey has shown us that non-symmetric systems are the mathematical signature of directionality and irreversibility in the universe. They appear when we model the flow of air and water, the response of materials to complex forces, and the webs of influence that structure our world. The solvers we have studied are therefore essential tools for a vast range of modern science and engineering. And their reach extends further still: by simply replacing real dot products with Hermitian inner products, these same algorithms can be used to solve problems in the complex domain, opening the door to modeling wave phenomena in quantum mechanics and electromagnetics. Far from being a mathematical curiosity, non-symmetry is a fundamental feature of the world we seek to understand and engineer.