Implicitly Restarted Arnoldi Method (IRAM)


The repeated application of (10) for every node, and the application of a similar equation for hy leads to an eigenvector problem of a large system of linear equations. The Finite Difference Mode Solver uses the Implicitly Restarted Arnoldi Method as described in Ref. [2] to find the eigenvectors of this system, and thereby find the modes of the waveguide.
The Implicitly Restarted Arnoldi Method looks for the modes inside a Krylov Subspace. This subspace is constructed from the mode operator, and from an arbitrary (could be random) starting magnetic field, here described as a vector, v1 . The Krylov subspace

Km ( Q, v1 )is the space spanned by powers of the Q when applied to the initial vector v1 :

BPM - Equation 12

The IRAM constructs a basis for the m + 1 dimensional Krylov subspace by the Gram-Schmidt process. For example, a basis for a 2 dimensional Krylov subspace, K1 ( Q, v1 ) , can be constructed from v1 and Q v1 . A vector perpendicular to  vcan be found by subtracting the v1   component from Qv1  :

BPM - Equation 13

This vector is parallel to the second basis vector, v2 , which is found once
w2 is divided by its length, BPM - Equation a . It is easy to show that this length is

BPM - Equation 14

The above inner product v of with Q occurs frequently, so the matrix components
are defined

BPM - Equation 15

With this definition, the effect of including another power of Q can be written in terms
of the new basis vectors

BPM - Equation 16

Following the Gram-Schmidt process, a vector in the Krylov subspace that is perpendicular to both v1 and v2 can be constructed from subtracting the projections of Qvonto the first two basis vectors

BPM - Equation 17

Normalizing w3 as before, and rearranging terms shows how the space generated

from Q2 vcan be represented in terms of the basis vectors

BPM - Equation 18

The Gram-Schmidt process can be followed in as many steps as necessary to create a basis for the Krylov subspace.

An important point to note is that the dimension of Km ( Q, v1 ) , (and therefore the number of basis vectors needed to span it), might be less than m + 1. Let’s suppose that the mode solving problem was asked to find the first k modes, where k < m. Let us also suppose that the initial vector v1 happened to be a linear combination of these k modes. Modes are solutions to the eigenproblem of equation (11). Let the ith mode field be represented by ei  , and let its propagation constant (eigenvalue) be βi2  . Call the subspace spanned by those k modes Sk , to represent a subspace of k dimensions. Then

BPM - Equation 19

When the second vector to span the Krylov subspace is calculated, it is apparent it is

BPM - Equation 20

also a linear combination of the k eigenvectors, and is therefore still in the subspace, Sk . It is easy to see that this situation is not changed by multiplying by Q again, and so

BPM - Equation 21

for any power of Q . (j  can be 1, 2, … up to ∞ ). Therefore the dimension of this Krylov subspace, Km ( Q, v1 ) is fewer than the m + 1 vectors in the span (12). For this choice of v1 , the Krylov subspace has dimension k ( < m), and some of the vectors in (12) have linear dependence on others.

When vis a member of the subspace Sk  , the Gram-Schimdt process must terminate at k, since having m basis vectors is not possible in Sk  . The termination will be seen at the kth step of the process, in trying to calculate a vector perpendicular to the previous k basis vectors. Equation (17) showed the generation of the 3rd vector, the kth will be formed from

BPM - Equation 22

When this is done with v1 ∈ Sk , the wk + 1 must be zero, and the Gram-Schmidt process will terminate. Because of the termination, the search for modes will take place in a space of dimension k ( 1, 2, maybe 10) rather than in the space of the vectors themselves. The vectors have dimension equal to the number of nodes in the mesh of Fig. 1. This is typically some number on the order of 105 or 106. This termination greatly simplifies the problem!

Of course, a random choice of v1 is unlikely to be a member of Sk , but the IRAM strategy is to use iterative processes that effectively change the initial choice of vto a vector inside (or not far from) the subspace Sk .

In the FD tab of the Solver Parameters dialog box of OptiMode, both k and m are defined as shown in Fig. 2.

BPM - Figure 2 Solver Parameter dialog box of OptiMode displaying the tab for the Finite Difference mode solver (FD). Dimension of Krylov Subspace (m) and number of modes to find (k) are defined here. In the Auto setting, m = 2k with m having a minimum of 4. m can be set by the user by unselecting the Auto check box.

Figure 2: Solver Parameter dialog box of OptiMode displaying the tab for the Finite Difference mode solver (FD). Dimension of Krylov Subspace (m) and number of modes to find (k) are defined here. In the Auto setting, m = 2k with m having a minimum of 4. m can be set by the user by unselecting the Auto check box.

Note that when the Dimension of Krylov Subspace is in its Auto setting, the dimension of the Krylov subspace (m) will be made double the dimension of the number of modes desired (k). In this example, the IRAM will construct a basis for a Krylov subspace of dimension m. This will involve making m copies of the field initially found on the defined mesh. If this mesh is already very large, there is a possibility to create an out-of-memory error, because the mode solver requires at least m times the storage required for a single mode.

The first choice of v1   will probably not be anywhere near Sk  , so the IRAM begins a process of shifted Q-R iteration to move away from unwanted eigenvalues and closer to the desired ones. When the Q-R iteration process is finished, The algorithm restarts (the R in IRAM). It reduces the dimension of the working space to k, but the new subspace is different because the QR iterations indirectly change the starting vector, v1 . The result is usually movement toward the desired solution space Sk  .

This leaves a new set of vectors spanning a new k dimensional space, and the algorithm starts the Arnoldi process again to build the dimension up to m. Then another Q-R process will start to move closer to Sk  . See Ref [2] for details. If you watch the computer’s Task Manager while the algorithm runs, you might see the used RAM cycle as the dimension of the space cycles from k to m and back. After many iterations the stopping condition will be met, and at this point the generalization of equation (18) will appear as

BPM - Equation 23

If the v1   is now inside Sk  , the last coefficient in the sum will be zero. This is evident from the definition of hij . If v1   really is inside Sk  , Qvk is also inside Sk   and it must be perpendicular to vk + 1 , making hk + 1, K = 0 . In the practical implementation of the mode solving problem, hk + 1, k won’t be exactly zero, but the IRAM process finishes when it is small enough to be negligible. Once the problem is in the form of (23) with the last term very small, it is easy to find the modes and propagation constants.

To see this, make the following change in notation. Define an n x k matrix columns of which are made from the k Arnoldi basis vectors of size n x 1, Vk  , the vk

BPM - Equation 24

and form a k x k matrix from the hij  coefficients of equation (15). Now consider QVk . Equations (16) and (18) show how the first two columns of QVk  are formed in terms of hij  and the vk   found in the first three columns of Vk  . The reader can find the following columns by following the Gram-Schmidt process further. Note that Hk  has all zero elements for those below the first subdiagonal, because the ith column of QVk is a linear combination of the first i + 1 basis vectors vk . (A matrix of this form is called “Hessian”).

BPM - Equation 25

All of the columns of QVk , except for the last one, are calculated by Vk Hk . The last column (kth column), according to (23) will need to include a term vk + 1 . This isn’t available in Vk , so the relation is made complete by adding the last column explicitly:

BPM - Equation 26

ek is the kth unit column vector in a k dimensional space, and the transpose is ekT = ( 0, 0, 0, …, 1 ). Putting this row vector after the column vector vk + 1 serves to put vk + 1 in the last column of a k x k matrix with zeros elsewhere. The second term on the right adds the missing column in the relation. From the above discussion, the hk + 1, k is expected to be small, and it becomes a measure of the accuracy of the solution.

The k modes we are looking for are inside Sk . These modes can be found by finding the eigenvectors of the matrix Hk . This greatly simplifies the problem. The matrix Q is n x n, where n is the number of nodes in the finite difference mesh. n is typically

on the order of 105 or 106. On the other hand Hk is k x k, and k is probably of order of 10 or so. The reduction to (26) simplifies the problem by many orders of magnitude.

The k modes inside Sk are a linear combination of the basis for that space. Let a vector, of size k x 1, s , represent an arbitrary vector in this space

BPM - Equation 27

Now let s be an eigenvector of Hk , and let the corresponding eigenvalue be θ .

BPM - Equation 28

It is easy to show that, within an accuracy specified by hk + 1, k , the x
defined in (27) is an eigenvector of Q with eigenvalue θ :

BPM - Equation 29

If the last term in (23) is zero, then ( x, θ ) is an exact eigenpair solution for Q . If the last
term is not exactly zero, then the pair ( x, θ ) is called a Ritz Pair, and the norm of (29) is a
measure of how accurately the Ritz Pair approximates the true eigenvector.

BPM - Equation 30

The number the above equation (30) evaluates to is called the Ritz estimate of the
Ritz Pair, ( x, θ ) .

The IRAM will stop when the user specified tolerance has been met, or the maximum number of steps has been exceeded. The tolerance is calculated from the user specified number in the Tolerance field of the Solver Parameters dialog box shown in Fig. 2. The maximum number of steps is from the field called Max. Steps. Letting the tolerance be t , the IRAM will stop when the following condition is met:

BPM - Equation 31

where εM is the machine precision, or when Max. Steps is exceeded.

When the mode solver is running, while the iterations are taking place, the value in (31) is
exceeded, but it is still informative to observe at least one of the Ritz vectors. In a long
calculation, the user can see if the mode solver is approaching a mode solution or not. If not, the simulation can be stopped, and corrective action taken.

OptiMode calculates one of the Ritz vectors by (27) while the IRAM is going through its iterations. This is the field displayed in the simulator while the calculation is taking place. You can select which Ritz vector to view with the entry in the field Preview Mode field in the Solver Parameters dialog box (Fig. 2). Selecting 0 in this field will calculate and display the first Ritz vector, the current approximation to the fundamental mode. Entering 1 in this field will calculate and display the second Ritz vector, an approximation to the next higher mode, and so on. The fields don’t change very much from one iteration to the next, so to save execution time you might not want to calculate and view a Ritz vector for every iteration. If you want to observe every 10th iteration, enter 10 in the Skip Value field.