Quasi-Newton Methods for Partitioned Simulation of Fluid–Structure Interaction Reviewed in the Generalized Broyden Framework

Fluid–structure interaction simulations can be performed in a partitioned way, by coupling a flow solver with a structural solver. However, Gauss–Seidel iterations between these solvers without additional stabilization efforts will converge slowly or not at all under common conditions such as an incompressible fluid and a high added mass. Quasi-Newton methods can then stabilize and accelerate the coupling iterations, while still using the solvers as black boxes and only accessing data at the fluid–structure interface. In this review, the IQN-ILS, IQN-MVJ, IBQN-LS, MVQN, IQN-IMVLS and IQN-ILSM methods are reformulated in the generalized Broyden framework to illustrate their similarities and differences. Also related coupling techniques are reviewed and a performance comparison is provided where available.


Introduction
Fluid-Structure Interaction (FSI) is the interaction of a fluid with a moving or deforming structure and occurs in many different branches of engineering. In mechanical engineering, the blades of a wind turbine deform due to their interaction with the wind [1,2]. Also Flow-Induced Vibration (FIV) can occur, for example in tube bundles with external flow, leading to leakage or even rupture of the tubes [3]. In civil engineering, there are interactions between wind flow and bridges [4], silos [5], tents [6] and many other structures. In the biomedical field, heart valves and arteries are flexible structures that interact with the blood flow [7,8].
As fluid-structure interaction is a multi-physics problem, complex phenomena can occur and numerical simulations are frequently used for the analysis. These numerical simulations of FSI can be performed in a monolithic or partitioned way. Monolithic codes solve the equations for the fluid and for the structure simultaneously, for example with a Newton-Raphson procedure [9] or multigrid method [10]. In this review, the focus is on the partitioned approach, as it allows to reuse mature and optimized codes to solve the subproblems.
Among the partitioned approaches, one can distinguish between the weakly and strongly coupled techniques. Weakly coupled techniques (also called explicit or loose coupling) solve the flow equations and the structure equations only once per time step [11,12]. Consequently, the equilibrium between the fluid and the structure is not satisfied exactly. In this context, the term equilibrium refers to the equality of velocities and forces on the fluid-structure interface. These weakly coupled techniques are typically suitable for aeroelastic simulations with light and compressible fluids [13], but specific schemes can also be applied with dense and incompressible fluids [14].
By contrast, strongly coupled techniques (or implicit coupling) use coupling iterations between the flow solver and structure solver to enforce the equilibrium at the fluid-structure interface up to a convergence tolerance in a steady simulation or in each time step of an unsteady one [15][16][17]. As a result, the flow problem and structure problem are solved multiple times per (time) step. Obviously, these coupling iterations increase the computational cost, but the cost per coupling iteration normally decreases during the iterations within a (time) step as the change per iteration decreases. In the remainder of the paper an unsteady simulation will be assumed; a steady simulation is then a special case with one large time step.
An important parameter for the choice between weakly and strongly coupled techniques, but also for the stability of the coupling iterations in several strongly coupled techniques, is the ratio of the added mass to the mass of the structure [18]. The added mass is the apparent increase in mass of the structure due to the fluid that is displaced by the motion of this structure. Physically, it is influenced by the shape of the fluid domain and the density of the fluid. Numerically, also the time step size determines its effect, but this effect depends on whether the fluid is compressible or not. FSI problems with a compressible fluid can always be stabilized as long as the time step is sufficiently small, regardless of the ratio of the apparent added mass to the structural mass. However, for an incompressible fluid, stability cannot be obtained by decreasing the time step size [13,15]. On the contrary, for incompressible flows in combination with flexible structures, decreasing the time step size may even increase the instability [15,[19][20][21]. For example, for the simulation of an elastic panel clamped at both ends and adjacent to a semi-infinite fluid domain, the added mass effect of a compressible fluid is proportional to the time step size, while for incompressible fluids it approaches a constant as the time step size decreases [13].
Especially for an incompressible fluid, many cases have a high added mass, e.g., blood flow in a vascular system [22], vibrations of tube bundles in lead-bismuth eutectic [23] or flutter of a slender cylinder in axial flow [24]. For these cases, the straightforward iteration between flow and structure solver within a time step will typically converge very slowly, or not at all, if no additional stabilization efforts are implemented. In this work, the focus is on techniques which consider the solvers as black boxes, as this is typically the case in a partitioned approach. Then, stabilization methods that alter one of the solvers, e.g., including an approximate added mass operator in the structural solver as in [25], are not possible. To stabilize and accelerate the convergence of coupling iterations with black box solvers, quasi-Newton methods have been developed in the fluid-structure interaction community. These methods will be reviewed in this work using the generalized Broyden framework.
The remainder of this review paper is structured as follows. First the FSI problem is posed and the necessary notation is introduced in Sect. 2. Then, the most basic solution approach is discussed in Sect. 3, with focus on its shortcomings in terms of stability and convergence speed, and how they can be overcome by introducing Jacobian information. In Sect. 4, a general method to obtain these Jacobians is discussed, called generalized Broyden. Thereafter, in Sect. 5, different quasi-Newton techniques are discussed in detail, including the Interface Quasi-Newton technique with an approximation for the Inverse of the Jacobian from a Least-Squares model (IQN-ILS), the Interface Quasi-Newton technique with Multi-Vector Jacobian (IQN-MVJ), the Interface Block Quasi-Newton technique with approximations from Least-Squares models (IBQN-LS), the Multi-Vector update Quasi-Newton technique (MVQN), the Interface Quasi-Newton Implicit Multi-Vector Least-Squares (IQN-IMVLS) and the Interface Quasi-Newton algorithm with an approximation for the Inverse of the Jacobian from a Least-Squares model and additional Surrogate Model (IQN-ILSM). This section ends with further notes and extensions on these methods. Finally, some numerical results to compare the different techniques are provided in Sect. 6, followed by the conclusions in Sect. 7.

Formulation of the FSI Problem
An abstract fluid-structure interaction problem, as shown in Fig. 1, consists of the subdomains Ω f and Ω s , with the subscripts f and s denoting fluid and structure, respectively. The boundaries of the subdomains are denoted as Γ f = Ω f and Γ s = Ω s and the fluid-structure interface Γ i = Γ f ∩ Γ s is their common boundary.
Besides having to satisfy the flow and structure equations in the respective subdomains while taking into account the appropriate boundary conditions on Γ f ⧵ Γ i and on Γ s ⧵ Γ i , the solution of the FSI problem is also required to fulfill the equilibrium conditions on the fluid-structure interface Γ i . The equilibrium conditions on a no-slip fluid-structure interface are twofold. First, the equality of fluid and solid velocity on Γ i is needed (kinematic condition) where ⃗ v is the velocity vector in the fluid domain and ⃗ u the displacement vector in the structure domain. Remark that this equality also implies equal accelerations on the interface. Second, equal magnitude but opposite sense of traction on Γ i is required (dynamic condition) Fig. 1 The fluid subdomain Ω f , the structure subdomain Ω s , their boundaries Γ f and Γ s and the fluid-structure interface Γ i where ̄f ,s is the stress tensor in Ω f ,s and ⃗ n f ,s the unit normal vector that points outwards from the corresponding subdomain.
As this work discusses coupling techniques that consider the solvers as black boxes, only the variables on the fluid-structure interface Γ i are of interest. However, the discretization of this interface is often different in the flow and structure subdomains. Given the focus of this review on coupling techniques, it is assumed that an interpolation layer is wrapped around or included in one (or both) of the solvers, invisible to the implementation of the coupling technique. As a consequence, the discretized displacement on either side of the fluid-structure interface can be represented as a column array x ∈ ℝ n x ×1 containing all components of the displacement vector ⃗ u in each of the n p grid points on the interface.
with the first subscript referring to the grid point (1 to n p ) and the second one to the component (1 to d, with d the dimension).
Similarly, the pressure p and all components of the viscous traction vector ⃗ t in each load point (1 to n l ) on either side of the fluid-structure interface are grouped in a column array y ∈ ℝ n y ×1 also called load vector, with the same meaning of the subscripts as above. Note that the n l load points do not need to coincide with the discretization of the displacement. It is important that the pressure load p ⋅ ⃗ n and viscous traction ⃗ t are not added, but included individually into y , because the pressure is typically dominant and has to stay perpendicular to the surface, also when interpolation is performed. If pressure and viscous traction were added, the resulting interpolated vector would have a pressure contribution that is not necessarily perpendicular to the surface after interpolation, resulting in an artificial shear component that can be much larger than the physical shear component.
With the typical Dirichlet-Neumann decomposition of the FSI problem, the displacement (linked to the velocity through the time discretization in time-dependent problems) is imposed at the interface in the flow solver and a pressure and viscous traction distribution is applied on the interface in the structure solver. A flow solver with a deforming grid using the Arbitrary Lagrangian-Eulerian (ALE) frame of reference will be assumed for the explanation, but this can be replaced by other techniques, for example the combination of the ALE approach and the Chimera technique [26] to handle large body motions or non-conforming alternatives, such as Immersed Boundary Methods (IBM) [27] and Embedded Boundary Methods (EBM) [28], which can handle large deformations and even topology changes. The flow calculation in a coupling iteration within a time step can be written as This notation concisely represents several operations and hides the dependence on previous time steps and the variables in the fluid domain next to the interface, while emphasizing the dependence on the discretized displacement x of the fluid-structure interface. It represents the following actions. First, the discretized displacement is given to the flow solver and the fluid domain adjacent to the interface is adapted accordingly. Then, the flow equations are solved in the entire fluid domain, resulting in a new load distribution y on the interface. Similarly, the calculation of the structure is represented by the function As before, this expression hides the dependence on both the previous time steps and the variables in the structure domain next to the interface. It indicates that the fluid pressure and viscous traction distribution on the interface y is given to the structure code. Subsequently, that code calculates the displacement of the entire structure and thus also the new displacement x of the fluid-structure interface.
With these notations, the FSI problem is formulated as the system that has to be solved for x and y . This problem can be rewritten as the root-finding problem with unknowns x and y. Moreover, the system in Eq. (7) can be reduced by substituting one equation in the other. Commonly, the first line is substituted in the second, but the other way around is equally possible. In this way, the FSI problem is simplified to a smaller system of equations which has to be solved for x . The notation • refers to function composition, so S • F(x) is equivalent with S(F(x)) . This looks like a fixed-point equation for x , but can also be written as a root-finding problem with unknown x (5) y = F(x).
To write this more compactly, the residual operator R(⋅) is defined as with output r = R(x) . The FSI problem thus reduces to finding the x that fulfills In this section, we have presented two formulations of the FSI problem. The first is the complete system Eq. (7) with n x + n y unknowns. The second is the reduced system Eq. (9), which has the benefit of having only n x unknowns. This system has been written more compactly using the residual operator R resulting in Eq. (12). In the next sections, several methods are discussed to solve the FSI problem presented here, in one of both formulations.
Both the solver operators as well as the residual operator are typically nonlinear. Therefore, the FSI problem exhibits similarities with nonlinear root-finding problems. The main difference is that an FSI problem usually involves time stepping (except for steady cases), which means that a nonlinear system has to be solved in each time step. Therefore, within each time step, coupling iterations are performed until the solution is reached. The nonlinear systems in subsequent time steps are somehow related to each other, because the solver operators change only gradually in time. As the solution is typically continuous, the initial guess for x at the start of each time step can be obtained by extrapolating the solution from previous time steps [29].

Gauss-Seidel Scheme
In order to solve the FSI problem, Eq. (8) has to be solved in each time step. One of the basic methods to solve such a system of nonlinear equations is the block Gauss-Seidel scheme. In this block-iterate scheme, each of the nonlinear equations is solved for one of the unknowns consecutively, and each unknown is updated to its new value as soon as it becomes available.
Because, further on, it will become necessary to make a distinction between the output of one solver and the input of the next, a tilde symbol is introduced to indicate the output of a solver: Using the superscript k + 1 to indicate the current iteration, the block Gauss-Seidel scheme takes the following form (12) R(x) = 0.
The lastly calculated displacement vector x k is used as x k+1 , the input of the flow solver in the following iteration. Subsequently, this vector is used to calculate a new load vector ỹ k+1 which is thereafter used as input of the structure solver y k+1 , to calculate a new displacement vector x k+1 . This iteration scheme, in which the output of the flow and structure solver is passed unchanged to the structure and flow solver, respectively, is the most basic way to find an equilibrium and is also called Gauss-Seidel scheme or fixed point iteration scheme.
The final solution of the FSI problem Eq. (7) has to fulfill the kinematic Eq. (1) and dynamic equilibrium condition Eq. (2) up to a certain tolerance. This means that x and ỹ have to approach x and y , respectively, which is expressed by the convergence conditions Because the output of each of the solvers is passed unchanged to the other, this can also be written as which relates to the fixed point formulation in Eq. (9).
By eliminating the occurrence of the load vector y , the procedure can be simplified to Furthermore, with the use of the residual operator introduced in Eq. (11), the iteration scheme becomes which is considered converged once

Motivation for Using Quasi-Newton Methods
Unfortunately, the Gauss-Seidel scheme explained above is not unconditionally stable due to among others the added mass effect.
Many researchers have investigated the stability of Gauss-Seidel iterations. For example, its convergence behaviour has been studied based on a simple model problem with a single degree of freedom on the interface [30]. Some investigated the added-mass effect [31]. Many others have explored the case of blood flow through a simplified artery [15,32,33]. They observed that, besides an apparent fluid mass of the similar order of magnitude as the actual structural mass, also a decrease in the stiffness of the structure or increase in domain length has a destabilizing effect. A first attempt to mathematically analyze the stability was done through the determination of the maximum relaxation factor to obtain convergence [15,25].
Instead of looking at a single number, the stability of a Gauss-Seidel scheme for a simplified flexible tube model with Dirichlet-Neumann decomposition can be examined by splitting the error on the interface into Fourier modes [34]. The mentioned error is the difference between the correct interface displacement and the one in a Gauss-Seidel iteration, based on linearized equations and without taking the boundary conditions into account. In this way, the authors were able to identify which frequency components become unstable. The analysis was first performed for a tube wall without inertia [34] and thereafter repeated including inertia [35], which proved to stabilize the convergence behaviour.
From this analysis it can be deduced that only a limited number of modes of the interface displacement are unstable and that the lowest wave numbers have the highest amplification factor and are hence the most unstable ones. This observation is true for different combinations of parameter values. In other words, the divergence or slow convergence of Gauss-Seidel iterations are caused by a limited number of unstable and slowly converging modes corresponding to the lowest wave numbers. The physical explanation for this observation is shown in Fig. 2. The figure shows an axisymmetric tube, the wall of which is perturbed with two different wave numbers, while on the in-and outlet a zero pressure boundary condition is imposed. Initially, its cross section is constant and the incompressible fluid is at rest. In the upper part of Fig. 2, a low wave number perturbation is applied and, because the fluid is incompressible, it is accelerated globally resulting in large pressure variations. In the lower part, a higher wave number perturbation is applied and the fluid acceleration is confined to more local regions. As a consequence, the pressure variations are much smaller for higher wave numbers. The pressure variations in the lower part of Fig. 2 are even barely visible, because the same scale is used for both cases.
Although the above analysis was performed on a flexible tube, the results are more widely applicable to incompressible fluids with a fluid-structure density ratio around one. For example, [37] arrived at the same conclusions by examining the stability of Gauss-Seidel iterations for a semi-infinite open fluid domain bounded by a string or a beam.
In summary, Gauss-Seidel iterations are not suitable for incompressible fluid cases with high added mass, because there is a limited number of error modes that are unstable. In order to obtain a solution for these cases, the unstable modes have to be removed by another technique to (efficiently) achieve convergence. Based on the results from the Fourier decomposition, it follows that only the low Fig. 2 The pressure contours (in Pa) in an axisymmetric tube due to two displacements of the tube's wall with the same amplitude but a different wave number. Initially, the fluid is at rest and the tube has a constant cross-section and zero pressure at both ends. A displace-ment of the tube's wall with a low wave number (top) creates much larger pressure variations than a displacement with a high wave number (bottom). Only the difference between the two calculations and not the values as such are important [36] 1 3 wave number modes have to be stabilized, while the others can still be treated using Gauss-Seidel iteration. The next section explains that the stabilization of these modes is achieved by including derivative information, which is the basic principle behind quasi-Newton techniques.

Quasi-Newton Schemes
In order to overcome this limitation of the Gauss-Seidel iterations for problems with high added mass and an incompressible fluid, the quasi-Newton iteration scheme is adopted. To improve stability, one or both of the vectors x and ỹ are modified before passing them to the other solver. If only one solver output is adapted, it is usually the output of the structural solver S as in Fig. 3b, but the opposite is equally possible. Figure 3c shows a schematic representation of adapting both solver outputs.
In the remainder of this section we will first introduce the adaptation of one output where we will use the modification of the structural output as example. This scheme will be referred to as the residual formulation scheme. Thereafter, the adaptation of both solver outputs will be introduced. The corresponding scheme is called the block iteration scheme.

Residual Formulation Quasi-Newton Scheme
In this scheme, the output x k of the structural solver is modified to x k+1 which is subsequently used as input for the flow solver. The output ỹ k of the flow solver is passed unchanged to the structural solver. Therefore, the load vector y can be left out altogether, as was the case for Gauss-Seidel iterations. With the use of the residual operator, defined in Eq. (11), the residual r k+1 in iteration k + 1 is written as and as before convergence is reached when Eq. (19) is satisfied. The difference with the Gauss-Seidel scheme is that x k+1 is no longer equal to x k . The adaption of the displacement vector follows from the use of a Newton-Raphson approach to solve the root-finding problem Eq. (12). This method uses the Jacobian of the nonlinear equation, which (20) is denoted here by R ′ , to estimate the input x k+1 that will direct the residual to 0 by solving for x k+1 . Note that Gauss-Seidel iteration Eq. (18) is retrieved if R � (x k ) = −I . Likewise, relaxed Gauss-Seidel iteration is obtained if the Jacobian is − I.
Because both the flow and structure solvers are considered black box solvers, the Jacobians of F and S are not accessible and hence, neither is R ′ . Therefore, the Jacobian of the residual operator is approximated, resulting in a quasi-Newton method where R′ (x k ) is the approximated Jacobian.
As explained in the previous section, the instability of Gauss-Seidel iterations is caused by a limited set of modes, i.e., for the vectors x in a small subspace of ℝ n x ×1 . Consequently, an approximation of the complete Jacobian of the residual operator R is not required. An approximated Jacobian which takes care of these unstable modes and leaves the other modes unchanged is sufficient. Leaving some modes unchanged means that the quasi-Newton method will actually perform Gauss-Seidel iterations for those modes. Solving the linear system in the equation above can be avoided by approximating the inverse of the Jacobian directly and calculating the update of the displacement vector as To conclude this section, a new short-hand notation is introduced for the approximate (inverse) Jacobians. For a nonlinear function r = R(x) , the approximate Jacobian and inverse Jacobian are written as

Block Iteration Quasi-Newton Scheme
Instead of only adapting the output of one solver, it is also possible to adapt the output of both the flow and structure solver. Now both x k+1 and y k+1 are different from x k and ỹ k+1 , and, because the load vector is no longer passed unchanged, it is not possible to use the residual operator. The convergence conditions are again given by Eq. (15).
The modification of the output of the solvers is determined by applying block Newton-Raphson iterations to the rootfinding problem Eq. (8) with unknowns x and y where Δx and Δy are the updates for the input x and y of the flow and structure solvers, respectively. Further, I is the identity matrix and F ′ and S ′ are the Jacobians of the flow and structure equations. Note that the two identity matrices will have different dimensions if the size of x and y differ.
Starting from the displacement x k that was given as input to the flow solver in the previous coupling iteration, the displacement x k+1 = x k + Δx k is calculated by solving the system for Δx k .
Using the updated value x k+1 and after calling the flow solver to determine ỹ k+1 = F(x k+1 ) , the pressure and viscous traction distribution y k+1 = y k + Δy k is calculated by solving the analogous system for Δy k . Subsequently, the structure solver is called to determine x k+1 = S(y k+1 ).
Similar to the previous section, the Jacobians are not accessible, because the solvers are considered black boxes. Therefore, approximations denoted by F ′ and Ŝ′ are used instead. Note that here two normal Jacobians are required, one for each solver, whereas in the previous section only one inverse Jacobian was required, namely the inverse Jacobian of the residual operator.
Adopting the same short-hand for the approximated Jacobians as in the previous section results in the following notations

Approximating Jacobians
The previous section introduced quasi-Newton approaches to stabilize and at the same time accelerate the convergence of coupling iterations. These schemes adapt either one or both of the solver outputs before passing them on, resulting in respectively, a quasi-Newton system for the residual formulation of the FSI problem, or a block iteration quasi-Newton system. These systems each contain one or more approximate Jacobians. The residual formulation scheme requires an approximation for the inverse of the Jacobian R � (x k ) , which is denoted as k r x . The block iteration quasi-Newton scheme requires approximations for the Jacobians of the flow solver and structure solver, so k xỹ approximates F � (x k ) and k yx approximates S � (y k ) . In these notations the superscript k refers to the iteration in which the Jacobian has been approximated.
All of these approximate Jacobians can be created using the generalized Broyden method. In this section, we will explain this method for the construction of an approximate Jacobian of an arbitrary nonlinear function b = B(a) . For now, we leave out the added complexity of FSI problems, for which an approximate Jacobian has to be constructed in each time step. This will be explained in the next section. Here, we just have an iterative method, where in each iteration k the Jacobian The same technique can also be used to approximate the inverse Jacobian B � (a k ) −1 by k b a. Instead of immediately presenting the rather complex generalized Broyden equation, it is introduced step by step, in a way that better fits the quasi-Newton FSI explanations found in literature.

Satisfying the Secant Conditions
The core idea of any quasi-Newton Jacobian approximation is to use the nonlinear function input-output information from previous iterations. Indeed, an input a i resulting in a certain output b i is a piece of valuable information about the behavior of the black box function B , which can be used to approximate the Jacobian B � (a k ) by k a b . In the current iteration k + 1 , the inputs of all previous iterations are available, as well as the corresponding outputs The input-output info is stored and used in the form of differences between consecutive iterations, defined as The notation refers to the difference between previous iterations, in contrast to the Δ notation, which refers to the desired change or update that needs to be performed.
Each pair ( a i , b i ) is called the secant information at iterations i and is related to a secant line to the nonlinear function B . Therefore, it can be interpreted as a finite difference approximation for the Jacobian in the direction a i . Furthermore, each secant information pair has a corresponding secant equation: If the approximated Jacobian k a b meets this secant condition, it uses a finite difference approximation for the actual Jacobian in the direction of a i , with the input-output information of iterations i and i + 1 . This secant information is relevant only if the Jacobian stays more or less the same during the k iterations, which means that B(a) has to behave close to linearly in the neighbourhood of a k .
The idea is to construct k a b , so that it fulfills all the k secant equations. To write this compactly, the differences defined previously in Eq. (31) are stored in the matrices A k and B k as follows Now, the k secant conditions can be collected in the matrix equation With n a and n b being the length of the input and output vectors, this is a system of n b k scalar equations for n b n a unknowns (the elements of the matrix k b a ). The system is thus typically underdetermined ( k < n a ). In order to find a unique solution, the least-norm solution is sought, which is in this case defined as the smallest matrix in the Frobenius norm that satisfies all secant conditions. The solution is given as where A k + is the pseudo-inverse 1 (or Moore-Penrose inverse) of the rectangular matrix A k , defined as To calculate the pseudo-inverse, it is necessary that the columns of A k are linearly independent. For now, we will assume this is always the case and the issue of linear dependence of the secant information is addressed in detail in the discussion on filtering in Sect. 5.1.
The expression for the approximate Jacobian presented above is elegant and short, but not very intuitive. Therefore, a different approach to obtain the same expression is given below.
The purpose of the approximate Jacobian is to determine an estimated change in output Δb that corresponds to an arbitrary change in input Δa , by evaluating For this purpose, the secant information from the previous iterations is utilised in the following approach.
First, the arbitrary vector Δa is approximated as a linear combination of vectors a i , i.e.
with c ∈ ℝ k×1 a coefficient vector.
It follows from the secant information that an input difference a i corresponds to an output difference b i , for 0 ≤ i ≤ k − 1 . Therefore, and under the assumption that the linear behavior of B is locally dominant, it can be stated that a linear combination of vectors a i will correspond to the same linear combination of vectors b i . This principle allows to determine Δb as Finally, it remains to determine the coefficients c . The system in Eq. (38) is typically overdetermined. Hence, the leastsquares solution for c will be used, which can be obtained by solving the square system of normal equations Therefore, the coefficient vector is given as Using this to calculate Δb results in Comparison with Eq. (37) reveals the same Jacobian as determined before in Eq. (35).
Matrix-free implementation Some of the algorithms explained in Sect. 5 require the explicit construction of the Jacobian matrix, while for others only its product with a vector, e.g., k a b Δa , is required. This last set of algorithms allows matrix-free implementation, for which the Jacobian matrix Eq. (35) never has to be calculated explicitly in practice, nor is the explicit calculation of the pseudo-inverse defined in Eq. (36) needed. How this is achieved is explained here.
Equations (40) and (41) show that the product of the pseudo-inverse with a vector is in fact the solution of the normal equations, but solving the normal equations Eq. (40) becomes unstable if the number of columns in the matrix A k is rather high. A more robust method to calculate the pseudo-inverse uses the reduced or economy-size QR decomposition [38] of A k where Q k A ∈ ℝ n a ×k is a matrix with orthonormal columns and R k A ∈ ℝ k×k is an upper triangular matrix. 2 Applying this to the normal equations Eq. (40) and using the fact that the inverse of R k A exists because the columns of A k are linearly independent, results in Symbolically, this means that the pseudo-inverse can be written as R k A −1 Q k A T , but it should never be constructed or stored. Instead, the product of the pseudo-inverse with a vector can be calculated by first evaluating the right hand side of Eq. (44) and subsequently solving the system using back-substitution, as R k A is an upper triangular matrix. The complete procedure to efficiently determine Δb given Δa is summarized in Algorithm 1.
In the following, the notation with the pseudo-inverse will still be used. Nonetheless, it should be kept in mind that the actual calculation has to be done using QR decomposition and back-substitution, avoiding the calculation of the inverse of matrices as well as the construction of large dense square matrices.

Adding an Initial Estimate for the Jacobian
Assuming the columns A k are linearly independent, the above obtained approximated Jacobian B k A k + is of rank k.
For the current discussion, it is assumed that k ≪ n a . Therefore, the matrix is a low-rank Jacobian approximation, and has an image or range of dimension k and a nullspace of dimension n a − k . As a result, with regard to its product with an arbitrary Δa , only the part of Δa ∈ range(A k ) will have a non-zero result. This becomes clear from the definition of the pseudo inverse in Eq. (36). The part of Δa ⟂ range(A k ) falls in the nullspace of the approximated Jacobian and the product of this part with the Jacobian is therefore zero. In other words, the approximated Jacobian is zero in every direction that is not a linear combination of the directions a i , for 0 ≤ i ≤ k − 1 , encountered in the previous iterations. Nonetheless, a full rank approximation of the Jacobian may be required, e.g., when it is used in a quasi-Newton method according to the residual formulation scheme. If this is the case, the current Jacobian approximation, using the approximation based on secant conditions for the part of Δa ∈ range(A k ) , can be expanded with an initial estimate of the approximate Jacobian 0 a b for the remaining part Δa ⟂ range(A k ).
The splitting of Δa in these two parts is based on orthogonal projection and visualized in Fig. 4. The orthogonal projection of a vector Δa onto the range of range(A k ) is given by Fig. 4 The vector Δa is split into a part inside the range of A k and another part perpendicular to that range This is the part of Δa ∈ range(A k ) . Using the complementary projector or just calculating the difference of Δa and its orthogonal projection gives the part of Δa ⟂ range(A k ) . Refer to [39] for a more complete discussion of projectors. Moreover, note that, using the QR decomposition these two parts are given by Now, the Jacobian approximation based on secant information can be extended with an initial Jacobian 0 a b: The questions why and how an initial Jacobian can be added have been answered. What remains is the choice of its value. Often, the identity matrix is used, scaled with a factor, typically −1 or − , which corresponds to (relaxed) Gauss-Seidel iteration, as explained below Eq. (21). This is the simplest approach to obtaining a full rank Jacobian approximation and will also be used in Sect. 5. In the case a low rank approximation suffices, e.g., the Jacobians for block iteration quasi-Newton, 0 a b = 0 can be used, which means the second term disappears completely. In still other situations, a physics-based surrogate may be available to use as initial Jacobian. This approach may accelerate convergence, but is application-specific and will be discussed further in Sect. 5.6.

Generalized Broyden Method
Up to now, the approximation of the Jacobian k a b was determined such that it met all secant conditions. However, this is not the only way to use the secant information. Another often used method (although not in FSI) is to only require the approximated Jacobian to fulfill the latest secant equation. Therefore, the matrices A k and B k only contain the latest piece of secant information: For all vectors Δa ⟂ range(A k ) (i.e., Δa ⟂ a k−1 ), we want to use the previous Jacobian k−1 a b . In other words, the effect of the approximated Jacobian remains unchanged in all directions orthogonal to a k−1 . This is called the no-change condition, which can be written formally as To obtain a Jacobian approximation with these specifica- where A k and B k now contain only one column. This is a recursive expression for the approximated Jacobian. In fact, this is Broyden's original method, 3 to construct the approximate Jacobian [40]. It was developed in the sixties, to solve systems of nonlinear equations. Furthermore, Broyden's method can be generalized. Instead of using only one secant condition and the approximate Jacobian from the previous iteration, m secant conditions can be used in combination with the approximate Jacobian from m iterations ago. This gives rise to the generalized Broyden method with This equation for the approximate Jacobian in generalized Broyden, however, can also be obtained in a more formal way, namely as the unique matrix that satisfies a number of conditions. Two equivalent ways are described in [41].
First, the approximated Jacobian can be obtained as the only matrix that simultaneously satisfies the m secant conditions in Eq. (34) and the n a − m no-change conditions Secondly, it can be obtained as the unique matrix that satisfies the m secant conditions Eq. (34) and minimizes the difference with the approximate Jacobian from m iterations ago, i.e.
where the subscript F denotes the Frobenius norm.
Furthermore, previously discussed methods are retrieved by choosing certain values for the parameter m in the generalized Broyden method. For m = 1 , Broyden's original method is recovered, while for m = k the pure secant method from Sects. 4.1 and 4.2 is obtained.
The generalized Broyden method was established much later than Broyden's original method. The first extension to the original one in the eighties led to a rather complex modified Broyden method [42,43]. In the nineties, Eyert [44] simplified this method by removing some nonessential parameters, resulting in the generalized Broyden method presented here.
Around that same time, the connection between the generalized Broyden method and Anderson acceleration (or Anderson mixing) was discovered. Anderson acceleration [45] was introduced in the sixties to accelerate fixed-point iterations. Based on the work by Van Leuken [46], Eyert showed that Anderson acceleration is mathematically equivalent to generalized Broyden with m = k , i.e., the pure secant method introduced in Sects. 4.1 and 4.2. This is not immediately apparent due to the very different ideas on which Anderson and Broyden originally based their methods.
In partitioned FSI simulations, several variants of the generalized Broyden method are used to approximate Jacobians in quasi-Newton iterations. These techniques were developed independently from the older methods (Anderson, Broyden and generalized Broyden) and the correspondence to those methods was only discovered recently [47,48]. Because a nonlinear system of equations has to be solved in every time step of an FSI simulation, there are some particularities with respect to Jacobian approximation, such as the reuse of secant information from previous time steps, as well as the removal of old and irrelevant secant information. These topics are discussed in the next section.
Computational complexity and storage This section ends with a first look into the computational complexity to obtain and use these approximate Jacobians. For simplicity, it is assumed that a and b are both vectors of length n a . This is usually not true for the block methods, but n b is typically proportional to n a . Further, it is assumed that n a ≫ k , i.e., the length of the vectors are much larger than the number of secant pairs available. More details will be provided later on for the different FSI methods. No details about the number of operations will be given, only the complexity of the leadingorder term will be discussed.
At the basis of the generalized Broyden method is the economy-size QR decomposition of A k , which is used for determination of the pseudo-inverse of A k . This QR decomposition is typically done with Householder transformations, resulting in a complexity of O(n a m 2 ) , which is also the total complexity of the evaluation of the product of this pseudoinverse with a vector. Already, it can be noted that, in the case that m = k , the computational cost quickly rises relative to a low fixed value for m.
If the approximate Jacobian is only needed to calculate its product with a vector, its explicit construction can be avoided. In some algorithms of the next section, however, the approximate Jacobian is used explicitly. Then, the construction of this n a × n a matrix has a complexity of O(n 2 a m) . In addition, the n a × n a matrix requires a storage capacity O(n 2 a ) , which is a strong disadvantage of these select algorithms.
In other algorithms, it is possible to avoid this expensive construction and use a matrix-free method to multiply the approximate Jacobian with a vector, i.e., without large dense square matrices. In practice, this is done by evaluating the product using Eq. (47) and multiplying the factors within each term from right to left. Then the complexity of this evaluation is only O(n a m 2 ) , which is the complexity of performing the QR decomposition needed to evaluate the product of the pseudo-inverse of A k with a vector. Because only the secant information has to be stored, the storage requirements O(n a m) are lower as well.

Difference Between the Anderson and Broyden Approach
The previous part formulated the generalized Broyden method, in which the parameter m determines how the secant information from previous iterations is included. Setting m = k corresponds to the Anderson method, in which the approximated Jacobian is determined by imposing all secant equations directly. For m = 1 , Broyden's original method is retrieved, here simply referred to as the Broyden method, in which the approximated Jacobian only fulfills the latest secant equation and the secant information from the previous iterations is included indirectly by imposing no-change conditions. In this section, the difference in behaviour between these two extreme versions of the generalized Broyden method will be clarified. Consider for example the approximation of the Jacobian B � (a k ) by k a b , when previously three iterations have been performed ( k = 2 ). The matrices containing the differences between consecutive iterations are Without loss of generality, it is stated that where p and q are orthonormal vectors, and x and y real scalars.
In the Anderson approach ( m = k ), the approximated Jacobian is The QR decomposition of A 2 is given by With this decomposition, the pseudo-inverse is calculated Finally, the approximate Jacobian is given by In the Broyden approach ( m = 1 ), the approximated Jacobian is Note that the pseudo-inverse of a single column equals its transpose divided by its norm squared, such that The resulting approximated Jacobian is given by Comparing Eq. (61) with Eq. (64), it is clear that the Jacobian approximations are different. Their inequality is analyzed in Table 1 by looking at their product with particular vectors Δa. For a vector Δa equal to the lastly added difference a 1 , both approaches return the corresponding difference b 1 , as expected. If the one before last vector a 0 is supplied, the results are different. The Anderson method simply returns b 0 , as this method attempts to approximate Δa as closely as possible using the already available differences. In contrast, the Broyden method does not and returns a linear combination of b 1 and b 0 . This approach gives priority to the lastly determined difference a 1 and uses the corresponding b 1 for the orthogonal projection of Δa on that difference a 1 . For a result of the Broyden approach that lies along b 0 , a difference orthogonal to the last difference a 1 needs to be supplied. Finally, the result for a general vector is given, where u, v and w are arbitrary scalars and s a unit vector orthogonal to p and q.
Both methods decompose Δa in components along the previously determined vectors a i , 0 < i < k − 1 , and   multiply the respective components with the corresponding vectors b i . This is shown graphically in Fig. 5. The Anderson method projects Δa on all previously determined differences. Therefore, the remaining part is orthogonal to these differences. The Broyden method projects Δa first on the lastly obtained difference, and the leftover part on the one before last, and so on. Therefore, the remaining part is not necessarily orthogonal to these differences. It will, however, always be orthogonal to the last difference onto which the projection was made, i.e., the oldest difference.
The difference between the two methods is essential to how nonlinearities in the secant information are dealt with. In general B(a k ) is nonlinear and its Jacobian is not constant, therefore the secant information will also contain nonlinear effects, especially when the step a i is large. Because the Broyden method prioritizes more recent secant information, it effectively ignores these nonlinearities, while the Anderson method does not, as it wants to approximate Δa as closely as possible using all available differences. This can lead to instabilities in the Anderson method. However, the Broyden method will also neglect small linear information, slowing down the convergence speed. More details and a method to remove nonlinearities from the secant information to stabilize Anderson are found in [49].
In the FSI community, the Anderson method is referred to as the least-squares approach and the Broyden method can be linked to the multi-vector approach. However, in FSI, fulfilling only the most recent secant equation as in the original Broyden method is typically not done and the multi-vector algorithms fulfill the secant equations in the most recent time step, using no-change conditions for older time steps. So in fact, the multi-vector approach is a generalized Broyden method, as will be explained in the following section.

Quasi-Newton Methods for FSI
In Sect. 3, different quasi-Newton schemes have been introduced. They required approximate Jacobians, which could be determined in different ways using information from (65) I − a 0 a 0 + I − a 1 a 1 + Δa previous iterations, as explained in Sect. 4. Up to this point, the focus was on solving the nonlinear equations in each time step separately. From here on, the distinction between different time steps will be necessary. Therefore, the superscript n + 1 will be used to indicate the values from the current time step, meaning that these are the values that are currently calculated. This notation is similar to the superscript k + 1 , which indicates the current iteration. In this section, the IQN-ILS, IBQN-LS, IQN-MVJ, MVQN, IQN-IMVLS and IQN-ILSM techniques will be derived and analyzed in the generalized Broyden framework. These techniques for partitioned FSI simulation have several differences, as summarized in Table 2.
The first difference is whether they use only the interface displacement as variables (IQN-ILS, IQN-MVJ, IQN-IMVLS, IQN-ILSM) or whether they are block iteration quasi-Newton methods using both interface displacement and load (IBQN-LS, MVQN). In the former case, they solve Eq. (12), in the latter they use Eq. (8), as explained in Sect. 3.
The second difference is related to how time stepping is handled, as most FSI simulations are time-dependent to capture a vibration or other dynamic behaviour. Assuming the inputs and outputs of the q previous time steps are stored, one can either impose the secant conditions from all time steps (IQN-ILS, IBQN-LS) or only for the latest time step, combined with no-change conditions for previous time steps (IQN-MVJ, MVQN, IQN-IMVLS, IQN-ILSM). This second difference is thus related to the choice of the parameter m of generalized Broyden, as explained in Sect. 4. On the one hand, m can be set to ∞ (actually limited to q time steps), so all the info from q previous time steps is used together, without an old Jacobian, but only with an initial one to start the procedure. In fact, this corresponds with the Anderson approach and is typically termed least-squares approach in FSI. On the other hand, only the secant info from the current time step can be used, with an old approximate Jacobian that is the final one from the previous time step. This corresponds to m = k , called multi-vector, and is really generalized Broyden, and not one of the limiting cases (Anderson or Broyden). The third difference is the amount of memory required for the storage of the approximate Jacobian(s) and the computational time required for the calculations related to the quasi-Newton steps. This will be explained more in detail for each method below and will be summarized in Table 3.

IQN-ILS
IQN-ILS is the abbreviation for Interface Quasi-Newton technique with an approximation for the Inverse of the Jacobian from a Least-Squares model [50]. The IQN-ILS technique performs an update of the input for the flow solver in each coupling iteration, using an approximation for the inverse of the Jacobian of the residual operator, so and Δr k = 0 − r k . The approximation for the inverse Jacobian R �−1 (x k ) ≡ k r x can be obtained directly by following the method explained in Sect. 4.3 for k a b , with a = r and b = x . Because the approximation needs to be full rank for a working quasi-Newton method, an initial Jacobian k−m r x = −I is used, which is the Jacobian of Gauss-Seidel iteration, as explained below Eq. (21). In the literature, however, the approximation for the inverse of the Jacobian R �−1 (x k ) is usually rewritten using the identity r =x − x , giving where the operator k rx is constructed as explained in Sect. 4.3, with a = r and b =x . In this way of explaining, no initial Jacobian is used, so k−m rx = 0 . It is worth mentioning that if instead of the inverse, the Jacobian R � (x k ) is approximated, the Interface Quasi-Newton Least-Squares method (IQN-LS) is retrieved [51].
For an FSI simulation with a single step (e.g., a steady simulation), the generalized Broyden formula in Eq. (52) is used with m = k and without initial guess 0 rx = 0 , so with Note that k rx has at most rank k, while In a time-dependent simulation, the secant information from the q previous time steps can be reused. As notation, the previous time steps are indicated with n, n − 1 , … , n + 1 − q , for the time step that is being calculated the superscript n + 1 is omitted and only the superscript k is used. The matrices R [k] and X [k] are a concatenation of the matrices from the different time steps, giving In this way, the information from each time step is treated equally, except when linear dependencies occur, because these are then removed by filtering, as will be explained below. The method thus satisfies all available secant conditions, i.e., from time step n + 1 and the q previous time steps. Consequently, m is equal to the number of columns in R [k] and X [k] , which is if no filtering is applied. As no initial Jacobian k−m rx is used, the information from earlier time steps n < n + 1 − q is not considered. It is important to remark that the difference between the first r or x of a time step and the last one from the previous time step is not used. Only differences between vectors of the same time steps are taken into account. This approach ensures that the secant information matches with the meaning of the Jacobian that is being approximated, which is the derivative within a time step, and not between them.
The reuse parameter q has to be defined by the user. Reuse typically improves the performance, but too old data is no longer helping the convergence, and therefore an optimal value exists. In the literature, the existence of this parameter is often cited as a drawback of IQN-ILS, because the performance of the method would be sensitive to this parameter. However, by using filtering, the performance of the method is rendered rather insensitive to this parameter around the optimum, as is shown by numerical tests in Sect. 6 and in other work [52,53]. Moreover, the parameter q allows the user to control how many time steps can be considered relevant, which is important in cases with rapid changes from one time step to the next, e.g., with multi-phase flows [54].
Matrix-free implementation Equation (69) is a symbolic notation to write IQN-ILS in the generalized Broyden framework, but this matrix should never be constructed or stored in the computer's memory. One of the main benefits of IQN-ILS is its so-called matrix-free character, which means that no large square matrices need to be constructed or stored. The product of the approximation of the inverse of the Jacobian with Δr k = −r k in Eq. (67) is symbolically calculated as In practice, the product c k = R [k] + Δr k of the pseudo-inverse of R [k] with Δr k is calculated using the economy-size QR decomposition and back-substitution, resulting in as explained in Algorithm 1. The vector Δr k is thus written as a linear combination of the r i , resulting in coefficients c k . As each r i has a corresponding x i , the change in x corresponding with Δr k can be obtained by calculating X [k] c k . The complete procedure can be found in Algorithm 2, with a relaxation step with factor on line 7, for the case in which R [k] and X [k] do not have any columns, e.g., at the beginning of the simulation.
Using R [ Filtering When columns of R [k] are linearly dependent up to a tolerance f , the diagonal elements of R [k] R in Eq. (74) become small and this system can no longer be solved accurately. Hence, an essential component of IQN-ILS is filtering, especially when data from previous time steps is reused [50]. Columns of R [k] that are linearly dependent up to the tolerance f need to be removed together with the matching columns in X [k] . As the newest information is stored on the left-hand side in R [k] , a r i that is a linear combination of newer r j (j > i) is removed. Columns can be removed if

ii referring to a diagonal element of R [k]
R [55]. The advantage of the first approach is that the tolerance f can be set by perturbing x with smaller and smaller changes until the change in x is no longer smooth, but numerical noise. In this case, the tolerance f can be considered as a measure of how accurate the flow solver and structural solver are calculating their solution. This filtering procedure is shown step by step in Algorithm 3. Obviously, it will be difficult to obtain convergence of the coupling iterations to a level that is lower than f .
Alternative filtering approaches are algebraic QR filtering and POD filtering [55]. In the algebraic filtering method (QR2), a column is removed if the diagonal element R . In the POD filtering, the eigenvalues of the autocorrelation matrix of R [k] are used to truncate old data. The numerical tests in [55] showed that algebraic QR filtering worked better than POD filtering or filtering using was not performed and remains as an interesting future work. Because the latter is directly related to the solver tolerances themselves, as explained above, it has been chosen for this work. Another reason to do filtering is to limit the number of secant conditions in cases with few degrees of freedom on the interface. Typically, m ≪ n x , but with only few degrees of freedom on the interface, the oldest columns of R [k] and X [k] need to be removed such that there are at most n x columns, to avoid an overdetermined Jacobian.
Computational complexity and storage The additional storage required for the IQN-ILS method is the matrices R [k] and X [k] , both ∈ ℝ n x ×m . Temporary storage is necessary for Q [k] R ∈ ℝ n x ×m and R [k] R ∈ ℝ m×m , and the small vector c k ∈ ℝ m . The storage thus scales linearly with the number of degrees of freedom in the interface's discretization. Furthermore, m can be reduced compared to Eq. (72) due to filtering. A rule of thumb is that it is typically not beneficial to include more than 50 columns.
The economy-size QR decomposition of R [k] has at most a complexity of O(n x m 2 ) if the fast Givens method or the Householder method is used [38]. The matrix-vector product in the right-hand side of Eq. (74) has a computational complexity of O(n x m) and solving the triangular system a complexity of O(m 2 ) . Consequently, also the computational complexity scales linearly with n x and is limited. In numerical tests, the IQN-ILS algorithm normally accounts for less than 1% of the total CPU time.

IBQN-LS
IBQN-LS stands for Interface Block Quasi-Newton method with approximation of the Jacobians using Least-Squares models (initially called reduced-order models) [21]. It uses the formulation in Eq. (25) and solves Eq. (26) and Eq. (27) in turn for Δx k and Δy k . Low-rank approximations for F � (x k ) and S � (y k ) are constructed using the generalized Broyden method with m as in Eq. (72) and without initial value for the Jacobian, so like in IQN-ILS. The reuse of previous time steps and the filtering are also applied in the same way as explained above. Consequently, this technique enforces all secant conditions from the current and q previous time steps, for both the flow solver and the structural solver.
For the approximate Jacobian of the flow solver F� (x k ) , the generalized Broyden framework is applied using a = x and b =ỹ . For m as in Eq. (72) and without initial value for the Jacobian this can be written symbolically as with where the secant information from the q previous time steps is combined with that from the current time step A symbolic formulation of the approximation Ŝ′ in the generalized Broyden framework can be obtained in a similar way, using a = y and b =x.
A disadvantage of this technique is that two linear systems need to be solved in each coupling iteration. In the original version, the n x × n x and n y × n y matrices corresponding with these systems were explicitly constructed using the symbolic notations as in Eq. (76) and they were solved with a direct linear solver [21]. However, by adopting an iterative linear solver like GMRES, only a procedure to calculate the product of the approximate Jacobians with a vector is required [56]. In practice, the number of iterations for the iterative solver is close to the number of columns used for the approximate Jacobians. Alternatively, the Woodbury matrix identity can be used to obtain a closed expression for the update [57].
Matrix-free implementation This matrix-free procedure will be explained here for the flow solver. When the product of F� (x k ) ≡ k x y with a vector Δx k needs to be calculated during the iterative solution of Eqs. (26) or (27), this can symbolically be written as For the practical implementation, Algorithm 1 is followed and this computation is split in two parts by the introduction of a coefficient vector c k , giving with c k the solution of The last part is the least-squares solution to an overdetermined system that can be solved efficiently by calculating the economy-size QR decomposition, followed by using back-substitution.
To summarize this procedure, the Δx k is decomposed as a linear combination of the columns in X [k] , then the observation is made that columns in X [k] and Ỹ [k] with the same index form a secant pair, such that the result can be approximated as the same linear combination of the columns in Ỹ [k] , as shown in Eq. (80). The complete procedure can be found in Algorithm 4.
Computational complexity and storage Compared to IQN-ILS, IBQN-LS requires approximately twice the memory, as the data for two approximate Jacobians needs to be stored. Furthermore, even though the matrix-free procedure with the iterative linear solver is faster than explicit matrix construction and direct linear solver, the computing time is higher than for IQN-ILS, where none of this is required. Nevertheless, the time required for the coupling algorithm scales linearly with n x and n y and thus remains small compared to that of the actual solvers. The solution of the linear systems could be avoided by writing the solution to Eqs. (26) and (27) symbolically, using matrix inverses and applying the Woodbury matrix identity as in [57]. In this way, the size of the matrices that have to be inverted is reduced from n x and n y to m.

MVQN
MVQN is the abbreviation for Multi-Vector update Quasi-Newton [58]. This method is based on the IBQN-LS method and is even identical to it in the first time step. However, the differences appear when data from previous time steps is included. MVQN considers the current Jacobian as the sum of the Jacobian from the previous time step plus a rank-k update. This update is then determined by enforcing the secant conditions from the current time step and minimizing the Frobenius norm of the update. This coincides with the generalized Broyden method with m = k n+1 and the initial Jacobian equal to the one from the previous time step.
Considering again the approximate Jacobian of the flow solver F′ , the generalized Broyden framework is applied using a = x and b =ỹ . The value m is now k = k n+1 and the initial value for the Jacobian is the one for the previous time step, giving Using the definition of the pseudo-inverse in Eq. (36), this can be reformulated as which corresponds with the original formulation in [58]. However, analyzing this method is most straightforward when considering Eq. (81). If this approximate Jacobian is multiplied with a vector Δx , then the secant information from the most recent time step is used for the part for which it is available, i.e., within the column span of the matrix X k . The previous approximate Jacobian is only multiplied with the leftover part of Δx , i.e., the part orthogonal to the column span of X k . So, even if secant information from previous time steps is available for the first part of Δx , it will not be used. This relates to the difference between the leastsquares and multi-vector approach explained in Sect. 4.4. The matrix Ŝ′ is constructed in a similar way Computational complexity and storage The main benefits of this method are that no parameter q is required and that filtering with a tolerance f is typically less essential as only a relatively small number of secant conditions from the most recent time step are considered. However, this comes at a significant cost, as the matrices F′ and Ŝ′ are constructed and stored in memory, such that the algorithm scales with n 2 x and n 2 y , both in terms of memory use as in computational complexity. Combined with the linear systems that have to be solved in each coupling iteration, this becomes expensive compared to the actual solver for a reasonably large number of degrees of freedom on the interface (e.g., more than 10 4 ). A linearly scaling adaptation is the RandomiZed Multi-Vector Quasi-Newton method (MVQN-RZ) [57], which will be explained in more detail with the other linearly scaling multi-vector methods at the end of Sect. 5.4. [58].

IQN-MVJ
The IQN-MVJ method is an acronym for Interface Quasi-Newton with Multi-Vector Jacobian [59]. This method adopts the idea for reuse from previous time steps proposed in MVQN and transfers it from the block iteration to the residual formulation quasi-Newton scheme, i.e., with only the interface displacement as variable. It is thus a generalized Broyden method which satisfies the secant conditions from the current time step while using the Jacobian from the previous time step as initial value. IQN-MVJ is linked with IQN-ILS in the same way as MVQN is linked with IBQN-LS. Except for the explicit construction of the approximate Jacobian, this method is identical to IQN-ILS in the first time step.
In IQN-MVJ, the approximation for the inverse of the Jacobian 5 is thus constructed as with Using the definition of the pseudo-inverse in Eq. (36), this can be reformulated as which corresponds with the original formulation in [59].
Computational complexity and storage The main drawback of IQN-MVJ is that the square matrix with the approximation for the inverse of the Jacobian is constructed and stored in memory, such that computational cost and memory requirement scale with n 2 x . Consequently, like MVQN, this approach becomes expensive compared to the solvers for a reasonably large number of degrees of freedom on the interface (e.g., more than 10 4 ).
To avoid this scaling and achieve linear complexity in n x and at the same time attempt to avoid a reuse parameter q, a matrix-free version of IQN-MVJ has been developed, named IQN-MVJ-RS-SVD [60]. Thereto, the approximation in Eq. (86) is reformulated as To avoid storage of a square matrix, this can be written using a recursive formula starting from 0 rx = 0 at time 0. Obviously, this requires storage of matrices � X and R for each time step, which is beneficial as long as the total number of columns is significantly smaller than n x . However, after a high number of time steps, this becomes prohibitively expensive and therefore three strategies are proposed in which the simulation is split into so-called chunks consisting of q ′ time steps after which a restart is performed [60]. The authors concluded that the IQN-MVJ-RS-SVD algorithm with a Singular Value Decomposition (SVD) restart strategy is the most promising. The approximate inverse Jacobian at the end of a chunk is then truncated by performing an SVD and truncating the singular values below a tolerance ′ f . This truncated SVD is then the initial Jacobian for the following chunk. Furthermore, the SVD is efficiently updated at the end of each chunk. As opposed to the original IQN-MVJ, the IQN-MVJ-RS-SVD method has linear complexity in n x , but the implementation is more elaborate than for IQN-MVJ. Compared to IQN-ILS, which has a reuse parameter q and QR filter tolerance f and also has linear scaling in n x , this method now requires a chunk size parameter q ′ and SVD filter tolerance ′ f , but with the claim that the performance is less sensitive to these parameters.
Another multi-vector method that achieves linear scaling memory requirements is the algorithm MVQN-RZ [57], which employs the randomized SVD not only to avoid the explicit construction of large square matrices, but also to circumvent the recursive reconstruction of the interface Jacobian. This algorithm has been applied to both block iteration quasi-Newton techniques and residual formulation quasi-Newton techniques. In every coupling iteration, the complexity of MVQN-RZ is O(z 2 n x + k 2 n x ) , where z is the number of decomposition modes. This compares to IQN-MVJ-RS-SVD with a complexity of O(z 2 n x + k 2 n x + Mkn x + k 4 ) in every coupling iteration, with z referring to the number of eigenvalues left after truncation and M the simulation chunk size. Additionally, this method requires an SVD update after every M steps, which has a complexity of O(Mz 2 n x ) . In providing these complexities, it is assumed that n x ≫ k, z, M.
and then k rx is called the approximation of the inverse Jacobian, but this notation is not used here to avoid confusion.

IQN-IMVLS
To mitigate the quadratic scaling in n x of IQN-MVJ, the IQN-IMVLS (Interface Quasi-Newton Implicit Multi-Vector Least-Squares) method has been developed with linear complexity in n x [53]. The first observation is that the factor � X k − n rx R k in Eq. (86) can be updated by adding one additional column in each coupling iteration, instead of recomputing it entirely so only the product n rx r k−1 needs to be evaluated. Determining the product of n rx with a vector requires a procedure to calculate a matrix-vector product. It is proven in [53] that n rx in Eq. (86) can be reformulated as if the initial Jacobian is assumed to be zero. This recursive formulation can be truncated after q terms, giving In IQN-IMVLS the reuse parameter q typically does not exhibit an optimum in terms of performance. Instead, both performance and computational cost grow with increasing values of q as the method converges towards the IQN-MVJ approach.
Computational complexity and storage As can be observed in Eq. (92), this procedure requires the storage of R j T for the q previous time steps. If the inverse of R j T R j is calculated via the LU decomposition using partial pivoting with row interchanges [38], then this scales with n x , but has slightly less robustness to bad conditioning than the Householder QR approach. As a result, the complete procedure has linear complexity in n x , like IQN-ILS. In addition, the QR decomposition is only applied on the secant information from the most recent time step, as opposed to the matrix with the secant information from all time steps in IQN-ILS. As not all secant information is combined into one matrix, the sensitivity to (almost) linear dependencies is smaller. Furthermore, no restart is required, but the implementation is a bit more involved than IQN-ILS or IQN-MVJ. It was also observed in [53] that including the secant information from the previous time step in X and R as well can accelerate the convergence, especially at the beginning of a time step. The complete procedure can be found in Algorithm 7.

IQN-ILSM
All techniques mentioned above use the flow solver and structural solver as black boxes, while sometimes the user has additional insight into the behaviour of the problem. This additional information can be incorporated to accelerate the convergence of the coupling iterations using the Interface Quasi-Newton algorithm with an approximation for the Inverse of the Jacobian from a Least-Squares model and additional Surrogate Model method (IQN-ILSM) [61]. In this technique, the secant information is combined with a so-called surrogate model, which behaves similarly to the actual solvers but is significantly faster. The origins of this technique can be found in the FreQ-LeSS algorithm for freesurface calculation, where secant conditions from previous flow solver iterations was combined with an analytical model of the problem [62,63]. The surrogate model is denoted as R s and the subscript s will be used to denote quantities related to the surrogate model. The inverse Jacobian of R s with respect to x is referred to as the surrogate Jacobian, which is assumed to stay the same during the entire time step. A procedure to calculate the product of this matrix with a vector is sufficient, without requiring construction and storage. To emphasize this matrix-free aspect, the surrogate Jacobian is represented by a function rxs (⋅) . Furthermore, this surrogate Jacobian can be either full-rank or only a low-rank approximation, with column space R s . At the end of each time step, the surrogate model is synchronized with the original model by interpolating the solution from the latter to the former. This approach avoids large discrepancies between the original model and the surrogate model after some time steps.
To gain insight in the IQN-ILSM technique, Eq. (75b) should first be revisited. This equation shows how IQN-ILS splits Δr k in a part R k R k + Δr k in the column span of R k and a part (I − R k R k + )Δr k orthogonal to it, using secant conditions for the former and Gauss-Seidel iteration for the latter. The term (I − R k R k + )Δr k for which no secant information is known can now be split once more into a part R s R s + I − R k R k + Δr k for which the surrogate model has information and the remainder I − R s R s for which Gauss-Seidel iteration is the best option. Obviously, the latter is zero if the surrogate Jacobian is full-rank. The split of Δr k and the approximate Jacobian used for each part can be written as Note that this equation is written in terms of X k rather than X k and that r x s (⋅) is used rather than rxs (⋅) . Here r x s refers to rxs − R s R s + and not to the full-rank Jacobian rxs − I as was the case in Eq. (68). Equation (93) shows that the secant conditions have the highest priority, followed by the surrogate Jacobian and then Gauss-Seidel iteration. In the first coupling iteration, there are no secant conditions yet, so the (93) surrogate Jacobian plays an important role. The contribution of the secant conditions will then become more significant during the coupling iterations as the column span of R k gradually increases. The expression in Eq. (93) can be simplified by using Jacobians with respect to x instead, giving These expressions can be expanded by including multiple surrogate models in a similar way. Several types of surrogate models can be considered. A first type is a coarse grid version of the original problem, using the same solvers. When using such a surrogate model, it is important to note that the secant information obtained on the original grid and on the coarse grid is not combined. The surrogate Jacobian thus also uses the coarse grid and interpolation is performed when it needs to be multiplied with a vector from the original grid. A second option is to use solvers with simplified physics, such as solvers neglecting viscosity or nonlinearity. If the simplified physics solvers have a known Jacobian, e.g., because they are analytical functions, then the surrogate Jacobian will typically be constructed and stored but it will be full-rank and only a relatively small matrix.
Another option for the surrogate is reuse from previous time steps in a time-dependent simulation. As opposed to the previously mentioned surrogate types, this surrogate model does not require the solution of a separate problem and no synchronization at the end of each time step, so it is essentially free. The reuse of q previous time steps corresponds to q nested surrogate models, with decreasing importance when for older time steps, giving with X i and R i containing the secant information from time step i. This equation can be condensed to When calculating the product ∏ n+1 i+1 (I − R j R j + )Δr k for each of the terms in the summation, the value in the previous term is stored and updated with one factor for the next term. Therefore the summation is in fact done in reverse order ( n + 1 � → n + 1 − q ). The IQN-ILSM algorithm can be found in Algorithm 8 and an efficient calculation of Δx s = rxs (⋅) in the case of reuse from previous time steps in Algorithm 9.
Computational complexity and storage While in IQN-ILS the QR decomposition is applied on the secant information from all time steps combined in each coupling iteration, IQN-ILSM with reuse as surrogate only applies the QR decomposition on the data from the current time step during the coupling iterations. The QR decomposition of the data from previous time steps is stored and does not need to be updated. This difference is reflected in the computational cost, which scales as O(n x (qk) 2 ) for IQN-ILS and O(n xk 2 ) for IQN-ILSM, with k the number of coupling iterations averaged over the time steps. Moreover, as IQN-ILSM does not combine the secant conditions from all time steps in a single matrix, it typically does not require filtering, as opposed to IQN-ILS.
By comparing Eq. (96) with Eq. (92), it can be observed that the IQN-ILSM method with reuse as surrogate is identical to the IQN-IMVLS method, except for some implementation aspects mentioned in Sect. 5.5. IQN-ILSM can thus be considered as a generalization of IQN-IMVLS such that not only reuse of secant information from previous time steps, but also physics-based surrogate models can be used, although it was not developed as such. Furthermore, the IQN-ILSM method can be interpreted as part of a larger class of methods which combine datadriven relations with physics-based knowledge [64].
In addition to providing a surrogate Jacobian to accelerate the convergence of the coupling iterations, the surrogate can also provide an initial guess for the coupling iterations. This prediction can be found in line 5 of Algorithm 8. When reuse from previous time steps is the surrogate model, this corresponds to linear extrapolation.

Other Algorithms
Besides the already mentioned techniques, there are many other variants, a selection of which is touched upon below. 6 This section ends with an outlook to the application of quasi-Newton coupling techniques outside of the field of fluid-structure interaction.
Jacobi iteration and high-performance computing The algorithms mentioned above are all using a sequential solution of the flow problem and structural problem. So, they can be considered related to Gauss-Seidel iteration, as opposed to Jacobi iteration which is characterized by simultaneous solution of both problems. For several of the quasi-Newton coupling techniques mentioned in this review, similar techniques are available based on Jacobi iteration [65]. For linear systems, Gauss-Seidel iteration inherently converge faster than Jacobi iteration, and this typically also holds for nonlinear systems, such as the FSI problem [66]. However, the higher number of iterations in a Jacobi-based method may be compensated by its better parallel scalability, although this is not guaranteed. Therefore, the difference in parallel scalability between both iterations types will be explained in the following lines. For Gauss-Seidel iteration, the flow and structure solvers are run sequentially, so both solvers can use all available CPU cores. However, the work load and the parallel efficiency of both solvers is typically very different, and often much higher for the flow solver. If the structure solver does not scale up well to a large number of cores, this leads to a low parallel efficiency, or analogously a number of idle cores. This problem is overcome by Jacobi iteration where the flow and structure solvers are calculating simultaneously. The cores are distributed over the two solvers such that both solvers require the same amount of calculation time, i.e., they are perfectly balanced, but this load balancing may not be trivial. Furthermore, specific variants of quasi-Newton methods have been developed for High-Performance Computing (HPC) such as the Compact Interface Quasi-Newton method (CIQN) [67], which is a parallel adaptation of IQN-ILS focused on efficiently combining partitions to realize a scalable implementation.
Multi-solver variants Also Multi-Solver (MS) versions of both IQN-ILS and IBQN-LS have been developed [68]. Once an FSI simulation can no longer be accelerated by increasing the number of cores per solver, the multi-solver algorithms can be applied for an additional speed-up. These multi-solver algorithms reduce the calculation time by running multiple instances of the flow solver and structural solver, while keeping the number of cores per solver constant and running each instance on one or more cluster nodes. One instance of the flow solver and of the structural solver perform coupling iterations like in the normal IQN-ILS or IBQN-LS algorithm. However, data from previous time steps is not reused directly as explained in Sect. 5.1, because the relation between the columns of R n and X n is only approximate at t n+1 . The additional instances of the flow solver and structural solver first recalculate the data from the previous time steps at the current time level, before including that data in a least-squares model. The columns of the matrix X n contain specific combinations of the degrees of freedom on the interface that accelerated the convergence of the coupling iterations in the previous time step. Hence, it is expected that knowing the difference of the output at t n+1 due to the same difference of the input as used at time level t n will improve the least-squares model for the approximate Jacobian.
Multi-level variants Furthermore, there exist Multi-Level (ML) versions of IQN-ILS and IBQN-LS [69]. Those could be considered similar to IQN-ILSM with a coarse grid surrogate model. However, the multi-level algorithms have important disadvantages compared to IQN-ILSM and therefore the IQN-ILSM algorithm is recommended. First, the multi-level algorithms combine the secant conditions obtained on all grids, which gives them all the same priority. By contrast, IQN-ILSM gives the highest priority to the finest grid, with a diminishing contribution from the coarser grid(s) as the coupling iterations on the finest grid converge. Furthermore, the multilevel algorithms interpolate the secant conditions obtained on the coarser grid(s) to the finest grid level and store it at the highest resolution, while IQN-ILSM stores the secant conditions with the resolution at which they have been calculated.
Aitken relaxation In addition to the quasi-Newton algorithms, Aitken relaxation [70][71][72] is frequently used for partitioned FSI simulations. This technique uses a dynamically varying scalar relaxation factor k for the Gauss-Seidel   The next input for R is thus a linear combination of the last output and the previous input. Therefore, the update of the interface's displacement is in the direction of the residual vector, as opposed to the update of the IQN-ILS method in Eq. (73). The value of k is calculated recursively as which can be interpreted as the secant method for scalars directly applied to vectors and projected on r k − r k−1 [70]. By combining Eqs. (97) and (98), it can be seen that the update of the interface's displacement is given by for k > 0 , which is similar to Eq. (73). If the Jacobian were created explicitly in the IQN-ILS algorithm and if the matrices R k and X k were limited to their newest column, Eq. (73) would yield Note the different location of the transpose sign in Eqs. (99) and Eq. (100). They are thus not identical because the coefficient of −r k is a scalar in the first equation and a matrix in the second one. Consequently, Aitken relaxation is different from IQN-ILS, even when the latter is restricted to one column in the matrices R k and X k . While Aitken relaxation uses a single relaxation factor for all interface degrees of freedom, IQN-ILS assigns a different value to each one based on a combination of the previously determined modes, i.e., the columns of R k and X k . Prediction In all the quasi-Newton algorithms, the coupling iterations begin from x 0 , which is an extrapolation or prediction. This can be a constant, linear or quadratic extrapolation based on x n , x n−1 , … . While there is an effect of the order of this prediction on the convergence of the coupling iterations, it is case dependent and it is thus difficult to state whether higher order is always faster [29]. In the IQN-ILSM algorithm, the surrogate model can also be used for a prediction, but this is also not always faster than a linear extrapolation [61].
Quasi-Newton methods with Robin-Neumann decomposition Instead of the typical Dirichlet-Neumann decomposition of the FSI problem, a Robin-Neumann decomposition can be used as well. This decomposition modifies the boundary condition in the flow solver to also include the pressure and traction forces. The effect is the introduction of Interface Artificial Compressibility (IAC), which allows to solve FSI (100) problems with enclosed incompressible fluids. To combine the idea of IAC with quasi-Newton methods a pressure correction has to be introduced to obtain corresponding inputs and outputs for the quasi-Newton technique [73]. However, by switching the order of the solvers and performing the quasi-Newton update on the pressure and traction forces instead of the displacements, this pressure correction can be avoided [74]. Applications in other fields Finally, it should be remarked that the quasi-Newton coupling techniques can also be used for coupled problems other than FSI. The main requirement for good results is significant interaction between both subproblems such that an approximate Jacobian stabilizes and accelerates the convergence compared to Gauss-Seidel iterations. The interaction should however not be so strong that an exact Jacobian satisfying all possible secant conditions is required to achieve convergence.
An example of a suitable problem is Conjugate Heat Transfer (CHT), where the partitioned spatial regions are each modelled by independent heat transfer codes and linked by temperature and flux matching conditions at the common boundaries [75]. A furnace radiation model can be coupled with a melt crystal growth model to investigate growth processes [76]. Another example is soil-structure interaction problems. The effect of excavations on the frame of a building can be studied by coupling a model of the soil's behaviour with a code for nonlinear structural dynamics [77]. In this case, the interaction between the models occurs at a relatively small number of points. Finally, they can even be used to calculate the combustion in a fluidized bed reactor under pressure [78]. The calculation of the carbon and oxygen concentrations is performed separately from the calculation of the temperature field. In this last example, there is only one domain and data are exchanged throughout the domain, as opposed to the other examples, where the domains do not overlap and data are only exchanged at the common boundary of the subdomains. Therefore, the cost of the coupling is no longer negligible compared to the solution of the subproblems, which necessitates coupling techniques with low computational complexity and low memory requirements.

Numerical Tests
Several comparisons of quasi-Newton techniques can be found in the literature [52,56,59]. This section aims at providing an idea about the relative performance of the described techniques using a test case that can be reproduced in a straightforward way, but without claim of general applicability. It focuses on the scaling with increasing number of degrees of freedom on the interface.

Test Case
The six quasi-Newton techniques presented in Sect. 5 are now evaluated on a one-dimensional (1D) flexible tube, through which an incompressible fluid flows. This test case runs quickly and still features the destabilizing added mass phenomenon [35].
The straight tube has a length , a circular cross-section with nominal inner diameter r 0 and wall thickness h. Its geometry is sketched in Fig. 6.
The ratio of the fluid density f to the structure density s is close to one, resulting in a large added mass. The structure has a modulus of elasticity E and a Poisson's ratio . The values of these parameters are reported in Table 4.
On the inlet, a pressure pulse of 1333.2 Pa is applied for a duration of 0.003 s. Due to the flexible tube wall, the pulse travels at a finite velocity, despite the incompressible fluid. The pressure at the outlet of the tube is atmospheric pressure, i.e., 0 Pa. The total simulation time is 0.01 s, divided in 100 time steps.
The flow solver solves the nonlinear continuity and momentum equation in which viscosity and gravity are neglected. The deformation of the tube wall is calculated by the structure solver, which only considers radial displacements, so the length of the tube stays constant. The tube is discretized in n p equal intervals. For the details regarding the governing equations and applied discretizations, the reader is referred to [35].
The FSI problem is solved with the open-source code CoCoNuT. This code allows to couple different software packages, which are treated as black boxes. It has the advantage of being modular and flexible in combining interpolators, solvers and coupling algorithms. Moreover, due to its comprehensive structure and implementation in Python, the code can be adapted to the user's own requirements, if needed. The most recent code can be found on GitHub in the repository https:// github. com/ pyfsi/ cocon ut/. The version used in this paper as well as scripts and result data are available for download on the Zenodo platform [79]. Both a flow and structure solver were written for the 1D flexible tube case described above and this code is equally available online.

Results
As the solvers are typically very expensive in FSI calculations, the focus is on the number of solver executions per time step. The term iteration will be used to denote a subsequent flow and structure solver call. Table 5 shows the number of iterations for the different quasi-Newton methods discussed above, as well as for a fixed relaxation factor and Aitken relaxation. Time data is included for completeness. A relative convergence tolerance is used, so that a time step is considered converged when ‖ ‖ r k‖ ‖2 < 10 −6‖ ‖ r 0‖ ‖2 , with r 0 the residual at the start of the time step.
Dynamic Aitken relaxation improves the convergence greatly compared to simple relaxation, but does not perform as well as the quasi-Newton methods. The block iteration and residual formulation methods perform similarly, whereas the least-squares techniques turn out to be somewhat more efficient compared to the multi-vector techniques for this test case. IQN-ILS (q = 0) corresponds to the method without reuse from previous time steps, but the IQN-ILS method performs best for q equal to 10. This shows the importance of reusing information from previous time steps. Note that the average number of coupling iterations per time step is not so sensitive to the value of q around the optimum. The indications (reuse) and (coarse) for the IQN-ILSM algorithm, refer respectively to the use of the previous time steps, and the use of an identical problem with half as many discretization intervals (without reuse) as surrogate model.
Next, Fig. 7 shows the memory requirements, for different number of discretization intervals n p . The data depict the peak memory usage as determined using the Python module guppy3. 7 The reported values are the total memory use minus that of the solvers. The results clearly show quadratically scaling memory requirements for the multi-vector techniques with explicit Jacobian construction, in contrast to the linear scaling for the methods with matrix-free implementation. Starting from 10 4 degrees of freedom, the memory requirements of MVQN and IQN-MVJ become much higher than those of the other coupling techniques. Furthermore, the block iteration methods require more memory compared to their corresponding residual formulation techniques, as two Jacobian matrices are approximated. Finally, IQN-IMVLS and IQN-ILSM (reuse) require more memory compared to the other matrix-free methods, because they retain information form every time step ( q = 100 ), instead of only from the last few ( q = 10 ). However, this could be reduced by applying a truncation after q time steps in the recursive formulas.
Lastly, Fig. 8 reports the coupling time, for different number of discretization intervals n p . The coupling time is calculated as the run time, excluding initialization, minus the time spent in both solvers and the interpolation routine. It is important to note that the code generating these results is not optimized for speed, such that care has to be taken when interpreting the results. Nonetheless, some relevant conclusions can be made. For example, remark that the quadratically scaling methods IQN-MVJ and MVQN require considerably more time than their linearly scaling counterparts. The same can be concluded for the block iterations methods, which are implemented here with an iterative procedure to obtain the quasi-Newton updates. Using the Woodbury matrix identity instead is expected to reduce the computational cost, but not to the extent that they will become cheaper than residual formulation techniques. Finally, the coupling time for the linearly scaling residual formulation techniques are in the order of 10 s for the investigated range of discretization intervals. These values remain low compared to the typical cost of solving the subproblems. Therefore, the difference in computational cost of the coupling for these methods is only moderately important.

Conclusion
The IQN-ILS, IBQN-LS, MVQN, IQN-MVJ, IQN-IMVLS and IQN-ILSM can all be reformulated as quasi-Newton techniques using generalized Broyden techniques for the approximate (inverse of the) Jacobian. In a time-dependent simulation with reuse from previous time steps, IQN-ILS and IBQN-LS enforce the secant conditions from all time steps in the same way, while MVQN, IQN-MVJ, IQN-IMVLS and IQN-ILSM only enforce the secant conditions from the latest time step directly. The block iteration quasi-Newton techniques like IBQN-LS and MVQN achieve similar performance in terms of number of coupling iterations per time step compared to the equivalent residual formulation quasi-Newton techniques IQN-ILS and IQN-MVJ, but the block iteration quasi-Newton techniques require more memory, more computational time for the coupling and a longer source code, entailing more involved implementation. Hence, it can be stated that, in general terms, it is better to use a residual formulation quasi-Newton method instead of the block iteration quasi-Newton techniques. For real-scale applications with a reasonable number of degrees of freedom on the interface, linear scaling of the coupling technique is essential, so MVQN and IQN-MVJ can only be used for smaller problems. The IQN-ILS technique has a short implementation and linear scaling thanks to the matrix-free procedure, but it requires filtering when reusing secant conditions from previous time steps. Finally, IQN-ILSM can be seen as a generalization of IQN-IMVLS, such that not only reuse from previous time steps can be applied but also other prior knowledge can be included to accelerate the convergence.