Free download functions of matrices: theory and computation

The use of numerical methods continues to expand rapidly. At their heart lie matrix computations. Written in a clear, expository style, it allows students and professionals to build confidence in themselves by putting the theory behind matrix computations into practice instantly.

Algorithms that allow students to work examples and write. Van Loan, published by Unknown which was released on Get Matrix Computations Books now! Matrix algorithms are at the core of scientific computing and are indispensable tools in most applications in engineering. This book offers a comprehensive and up-to-date treatment of modern methods in matrix computation. It uses a unified approach to direct and iterative methods for linear systems, least squares and eigenvalue problems.

When reality is modeled by computation, matrices are often the connection between the continuous physical world and the finite algorithmic one. Usually, the more detailed the model, the bigger the matrix, the better the answer, however, efficiency demands that every possible advantage be exploited.

The articles in this volume are. Applies matrix techniques to the solution of linear systems of equations and eigenvalue problems. Algorithms and computer implementation are presented, and the treatment of sparsity in large order systems and accuracy control are discussed in the light of practical applications.

Revised and updated, the third edition of Golub and Van Loan's classic text in computer science provides essential information about the mathematical background and algorithmic skills required for the production of numerical software. This new edition includes thoroughly revised chapters on matrix multiplication problems and parallel matrix computations, expanded treatment. The topics include representations, fundamental analysis, transformations of matrices, matrix equation solutions as well as matrix functions.

Mathematical Theory of Computation. Theory of computation. Elements of computation theory. The Theory of Computation. Computation Theory and Logic. Recommend Documents. Higham University of Manchester Ma Axler K. Your name. Conversely, when f is analytic we can obtain 1. However, Theorem 1. As an application of Theorem 1.

Miscellany In this section we give a selection of miscellaneous results that either are needed elsewhere in the book or are of independent interest. Applying Theorem 1. We noted in Section 1.

The next result shows that all matrices that commute with A are polynomials in A precisely when A is nonderogatory—that is, when no eigenvalue appears in more than one Jordan block in the Jordan canonical form of A. This result is a consequence of Theorem 1. The next result shows that a family of pairwise commuting matrices can be simultaneously unitarily triangularized. This corollary will be used in Section Related to Theorem 1. In this representation H can be taken to be Hermitian positive definite.

Note that another way of expressing Theorem 1. Interest in computing matrix functions grew rather slowly following the advent of the digital computer. As the histogram on page indicates, the literature expanded rapidly starting in the s, and interest in the theory and computation of matrix functions shows no signs of abating, spurred by the growing number of applications.

The original version of MATLAB included the capability to evaluate the exponential, the logarithm, and several other matrix functions. The availability of matrix functions in MATLAB and it competitors has undoubtedly encouraged the use of succinct, matrix function-based solutions to problems in science and engineering.

Notes and References The theory of functions of a matrix is treated in a number of books, of which several are of particular note. Almost every textbook on numerical analysis contains a treatment of polynomial interpolation for distinct nodes, including the Lagrange form 1. Theorems 1. We derived Theorem 1. The special case of Theorem 1.

These are more complicated and less explicit than 1. Problems The only way to learn mathematics is to do mathematics. Let Jk be the Jordan block 1. Use Theorem 1. Deduce the corresponding result for column sums. Use 1. What are the interpolation conditions 1. Obtain an explicit formula for f A. Show how to obtain the formula 1.

Prove the formula 1. Use this formula to derive the Sherman—Morrison formula B. Deduce the Cayley— Hamilton theorem. Is the commutativity condition necessary? Hint: Problem 1. Verify the Cauchy integral formula 1. There are in fact only two square roots of any form, as shown by Theorem 1. Describe all such A. Extension of Theorem 1.

Show that there is a unique square root X of A that is a primary matrix function of A and whose nonzero eigenvalues lie in the open right half-plane. Show that if A is real then X is real. Hint: use Theorem 1. Show further that if A is nonderogatory then X is a polynomial in A. Assume also that B and f A are nonsingular. What additional hypotheses are required for this proof? Give another proof of Theorem 1. Show that Corollary 1. Can 1. This was Putnam Problem B4.

This paper gives Cayley considerable claim to the honor of introducing the modern concept of matrix, although the name is due to Sylvester There results the important theorem that the latent roots of any function of a matrix are respectively the same functions of the latent roots of the matrix itself. All of the definitions except those of Weyr and Cipolla are essentially equivalent. FERRAR, Finite Matrices If one had to identify the two most important topics in a mathematics degree programme, they would have to be calculus and matrix theory.

Noncommutative multiplication underlies the whole of quantum theory and is at the core of some of the most exciting current research in both mathematics and physics.

Chapter 2 Applications Functions of matrices play an important role in many applications. We describe some examples in this chapter. Rearranging or reformulating an expression may remove or modify the f A term, and it is always worth exploring these possibilities. Here are two examples in which computation of a matrix square root can be avoided. The Cholesky factorization is generally much less expensive to compute than the square root.

However, the matrix exponential is explicitly used in certain methods—in particular the exponential integrators described in the next subsection. Exponential integrators are a broad class of methods for 2. They have attracted renewed interest since the late 37 2. The method 2. Nuclear Magnetic Resonance Two-dimensional nuclear magnetic resonance NMR spectroscopy is a tool for determining the structure and dynamics of molecules in solution. This relation is used in both directions: in simulations and testing to compute M t given R, and in the inverse problem to determine R from observed intensities.

The latter problem can be solved with a matrix logarithm evaluation, but in practice not all the mij are known and estimation methods, typically based on least squares approximations and requiring matrix exponential evaluations, are used.

Consider a time-homogeneous continuous-time Markov process in which individuals move among n states. The row sums of P are all 1, so P is a stochastic matrix.

Now consider a discrete-time Markov process with transition probability matrix P in which the transition probabilities are independent of time. If such a Q exists it is called a generator and P is said to be embeddable.

Every primary logarithm is complex, since it cannot have complex conjugate eigenvalues. Indeed, the standard inverse scaling and squaring method for the principal logarithm of a matrix requires the computation of a matrix root, as explained in Section The latter authors, who are motivated by credit risk models, address the problems that the principal root and principal logarithm of P may have the wrong sign patterns; for example, the root may have negative elements, in which case it is not a transition matrix.

They show how to optimally adjust these matrices to achieve the required properties—a process they term regularization. Their preferred method for obtaining transition matrices for short times is to regularize the appropriate matrix root.

Questions of existence and uniqueness of stochastic roots of stochastic matrices arise see Problem 7. Therefore the matrix exponential and logarithm are needed for these conversions.

We turn now to algebraic equations arising in control theory and the role played by the matrix sign function. The matrix sign function was originally introduced by Roberts [] in as a tool for solving the Lyapunov equation and the algebraic Riccati equation. Hence we can apply the sign function to 2. The system is consistent, by construction, and by rewriting 2. To summarize: the Riccati equation 2.

From 2. Moreover, writing Z in 2. Theorem 2. See Problem 2. The number of eigenvalues in more complicated regions of the complex plane can be counted by suitable sequences of matrix sign evaluations. The use of the matrix sign function as described in this section is particularly attractive in the context of high-performance computing.

Orthogonalization and the Orthogonal Procrustes Problem In many applications a matrix that should be orthogonal turns out not to be because of errors of measurement, truncation, or rounding. The question then arises how best to orthogonalize the matrix. One approach is to apply Gram—Schmidt orthogonalization or, equivalently, to compute the orthogonal factor in a QR factorization.

Nevertheless, there are strong connections with the matrix square root and the matrix sign function and so we devote Chapter 8 to the polar decomposition.

Optimal orthogonalization is used in a number of applications. In these applications the matrices A and B represent sets of experimental data, or multivariate samples, and it is necessary to determine whether the sets are equivalent up to rotation. An important variation of 2. Many other variations of the orthogonal Procrustes problem exist, including those involving two-sided transformations, permutation transformations, and symmetric transformations, but the solutions have weaker connections with the polar decomposition and with matrix functions.

Theoretical Particle Physics Lattice quantum chromodynamics QCD is a research area of physics that has in recent years made extensive use of matrix functions. QCD is a physical theory that describes the strong interactions between quarks as the constituents of matter. Lattice QCD formulates the theory on a four dimensional space-time lattice, and its numerical simulations currently occupy large amounts of high-performance computer time.

The system 2. Hence a key step in the solution of 2. In view of the huge dimension of H, this product must be computed without forming the dense matrix sign H. A variety of methods have been proposed for this computation, including Krylov techniques and methods based on rational matrix sign function approximations.

This identity converts the problem into one of estimating the diagonal elements of log A. Quantum Monte Carlo QMC simulations with the Hubbard model of particle interactions rely on a number of fundamental linear algebra operations. Other Matrix Functions One matrix function generates the need for another.

Thus eA may generate the need for log A and cos A for sin A. The inverse scaling and squaring algorithm for the matrix logarithm requires matrix square roots see Section Nonlinear Matrix Equations One reason for developing a theory of matrix functions is to aid the solution of nonlinear matrix equations. Ideally, closed form solutions can be found in terms of an appropriate matrix function. The algebraic Riccati equation 2. Here, any square root of AB can be taken.

The formula 2. Another expression for X is as the 1,2 block of sign A0 B0 , which follows from the sign-based solution of 2. As this example illustrates, there may be several ways to express the solutions to a matrix equation in terms of matrix functions.

The uniqueness follows from writing 2. Note that this approach leads directly to 2. One way of solving 2. The standard approach is to reduce 2. Using the theory of matrix square roots the equation 2. By using the Jordan canonical form it is usually possible to determine and classify all solutions as we did in Section 1. Finally, we note that nonlinear matrix equations provide useful test problems for optimization and nonlinear least squares solvers, especially when a reference solution can be computed by matrix function techniques.

The geometric mean has the properties see Problem 2. Pseudospectra Pseudospectra are not so much an application of matrix functions as objects with intimate connections to them. Pseudospectra provide a means for judging the sensitivity of the eigenvalues of a matrix to perturbations in the matrix elements. For example, the 0. More generally, pseudospectra have much to say about the behaviour of a matrix. It includes algorithms for the exponential, logarithm, square root, and trigonometric functions, all based on algorithms for matrices.

Sensitivity Analysis Ideally a numerical algorithm returns not only an approximate solution but also an estimate or bound for the error in that solution. A separate question, usually easier to answer, is how sensitive is the solution of the problem to perturbations in the data. Sensitivity is determined by the derivative of the function that maps the input data to the solution. Other Applications Finally, we describe some more speculative or less well established applications.

The function f agrees with ez up to terms in z of second order, and the idea is that f may bring more favourable stability properties than a polynomial or rational stability function when the eigenvalues of the Jacobian matrix vary greatly in magnitude and location. Schmitt uses the unscaled Denman—Beavers iteration 6.

Various algorithms are available for the solution of such problems. For a given p, the matrix sector function 2. Figure 2. However, a good numerical method for computing the matrix sector function is currently lacking. It can be used to obtain invariant subspaces in an analogous way as for the matrix sign function.

Note that this is the same mathematical idea as in Section 2. Bregman Divergences Matrix nearness problems ask for the distance from a given matrix to the nearest matrix with a certain property, and for that nearest matrix.

The nearest unitary matrix problem mentioned in Section 2. Thus Bregman divergences provide another application of matrix functions. In order to approximate Y at non-mesh points t it is necessary to interpolate the data, which can be done by standard procedures. However, if Y t has structure the interpolated matrices may not possess that structure.

The idea is to choose f so that the required structure is enforced. The expression 2. However, in these early papers nonunitary transformations are used, in contrast to the approach in Section 2. Early references for 2. Proofs of 2. The geometric mean 2. The factorization 2. Problems Though mathematics is much easier to watch than do, it is a most unrewarding spectator sport. Derive the formula 2.

Reconcile the fact that the initial value problem 2. Prove Theorem 2. Prove the relations 2. Problems 53 Whenever there is too much talk of applications, one can rest assured that the theory has very few of them. Existing techniques were often quite unreliable. For this reason, implementations are provided using only matrix inversion, multiplication, and addition.

In fact this is the case, although the presence of W often goes unrecognized. Computations with exact data are subject to rounding errors, and the rounding errors in an algorithm can often be interpreted as being equivalent to perturbations in the data, through the process of backward error analysis.

Therefore whether the data are exact or inexact, it is important to understand the sensitivity of matrix functions to perturbations in the data. Sensitivity is measured by condition numbers.

Most of the results of Section 3. We will use the condition number 3. Some care is needed in interpreting 3. The bound 3. An example is given in Section 5.

Usually, it is the relative condition number that is of interest, but it is more convenient to state results for the absolute condition number. To obtain explicit expressions analogous to 3. If we need to show the dependence on f we will write Lf X, E. When we want to refer to the mapping at X and not its value in a particular direction we will write L X.

In view of 3. From 3. See the works cited in the Notes and References. Depending on the context, 3. Theorem 3. The following theorem will also be very useful. For the rest of this section D denotes an open subset of R or C. We actually need only a special case of 3. A Theorem 3. Therefore since 3. Moreover, G X, E is a linear function of E. This follows from 3. Under the conditions of Theorem 3.

This is useful both in theory and in practice. Pm Proof. Then see the more general Problem 3. From B. Now consider a general function f. Problem 3. Note that Theorem 3. Corollary 3. Using 3. Bounding the Condition Number Our aim in this section is to develop bounds for the condition number. The bound of the theorem is obtained by maximizing over all the eigenvalues, using 3.

By specializing to the Frobenius norm we can obtain an upper bound for the condition number. By Corollary 3. A normal matrix is diagonalizable by a unitary similarity. The nonnormality term can of course be arbitrarily large. Computing or Estimating the Condition Number The essential problem in computing or estimating the absolute or relative condition number of f at X is to compute or estimate kL X k.

For the Frobenius norm the link 3. Algorithm 3. For large n, Algorithm 3. No analogous equalities to 3. Lemma 3. Hence, using 3. One of the uses of the condition number is to estimate the error in a computed result produced by a backward stable method.

We have omitted the normalizations to avoid cluttering the algorithm. An alternative to the power method for estimating the largest eigenvalue of a Hermitian matrix is the Lanczos algorithm. He shows that the Lanczos approach generates estimates at least as good as those from Algorithm 3. Turning to the 1-norm, we need the following algorithm. Key properties of Algorithm 3. Thus the algorithm is a very reliable means of estimating kAk1. We can apply Algorithm 3.

Advantages of Algorithm 3. Algorithms 3. A rough guide to the choice of t can be developed by balancing the truncation and rounding errors. Relative errors in the Frobenius norm for the finite difference approximation 3. In practice, topt is usually fairly close to minimizing the overall error. Figure 3. Then the underlying unitary transformations, exploiting the fact that G method is used to evaluate f at the triangular matrix. Several independent E can be used in order to get a better estimate.

For example, with two E the estimate is within a factor 5 of kL X kF with probability at least 0. The identity 3. Their proofs work with K X , which in this case has an explicit representation in terms of Kronecker products; see Problem 3.

That Corollary 3. In these analyses the assumption of a dominant eigenvalue is needed to guarantee convergence; for Algorithm 3. The power method in Algorithm 3. Prove 3.

In this chapter we survey a variety of techniques applicable to general functions f. Then we turn to methods based on similarity transformations, concentrating principally on the use of the Schur decomposition and evaluation of a function of a triangular matrix. Matrix iterations are an important tool for certain functions, such as matrix roots. We describe the cost of an algorithm in one of two ways. If, as is often the case, the algorithm is expressed at the matrix level, then we count the number of matrix operations.

For us, a method is usually a general technique that, when the details are worked out, can lead to more than one algorithm. Matrix Powers A basic requirement in some methods for matrix functions is to compute a power of a matrix. A sequence of successive powers A2 , A3 ,. Instead, repeated squaring can be used, as described in the following algorithm.

The initial while loop simply ensures that multiplication of a power of A with the identity matrix is avoided. Algorithm 4. This algorithm evaluates the polynomial 4.

In this case pm can be evaluated by explicitly forming each power of X. Note that while Algorithms 4. One drawback to Algorithm 4. In such situations the algorithm can be adapted in an obvious way to employ a factorization of pm into real linear and quadratic factors. The powers X 2 ,. Number of matrix multiplications required by the Paterson—Stockmeyer method and Algorithms 4.

Note that 4. As an 2 extreme case, this method evaluates Aq as Aq q , which clearly requires much less work than the previous two methods for large q, though for a single high power of A binary powering Algorithm 4. Table 4. For each m in the table it can be shown that both choices of s minimize 4.

The next theorem provides error bounds for three of the methods. Theorem 4. The computed polynomial pbm obtained by applying Algorithm 4. This is true even in the case of scalar x. Turning to Theorem 4. Error analysis for Algorithm 4. The main thing to note about the bound 4. We will not consider Algorithm 4. Taylor Series A basic tool for approximating matrix functions is the Taylor series. Truncation errors are bounded in the following result. Suppose f has the Taylor series expansion 4.

Moreover, the norm need not be consistent. For certain f this is straightforward. We illustrate using the cosine function. A2i , the bound of Theorem 4.

From Theorem 4. In fact, since 18! An advantage over polynomials is that they are better able than polynomials to mimic the behaviour of general nonlinear functions. However, a good approximation on the spectrum does not imply a good approximation overall. Therefore, to produce useful error bounds either rational approximations need to be applied to a restricted class of A or approximation errors must be analyzed at the matrix level. Let Rk,m denote the space of rational functions with numerator and denominator of degrees at most k and m, respectively.

These approximations are usually employed for Hermitian matrices only, so that error bounds for the scalar problem translate directly into error bounds at the matrix level. It is usually required that pkm and qkm have no common zeros, so that pkm and qkm are unique. The condition 4. Evaluating Rational Functions An important question is how to evaluate a given rational function at a matrix argument.

The obvious approach is to evaluate pkm and qkm by any of the methods discussed in Section 4. When the Paterson—Stockmeyer method is used to evaluate pmm A and qmm A , some savings can be made.

Referring to the description in Section 4. Consider now the continued fraction representation 4. This expansion can be evaluated at the matrix X in two ways: top down or bottom up.

This algorithm evaluates the continued fraction 4. Using bottom-up evaluation, rm X is evaluated as follows. The top-down evaluation is computationally expensive, but it is well suited to situations in which the aj and bj are independent of m and the whole sequence r1 X , r2 X ,. Another representation of rm in 4. The cost of evaluating 4. Since the stability depends very much on the function f , we delay further consideration until later sections on particular f.

Listing 4. Suppose the only error in the process is an error E in evaluating f B. This method may give inaccurate results if V is badly conditioned. For vector arguments, the function is applied to each component.

The diagonalization approach can be recommended only when the diagonalizing transformation is guaranteed to be well conditioned. This class includes orthogonal matrices, symmetric matrices, skew-symmetric matrices, and their complex analogues the unitary, Hermitian, and skew-Hermitian matrices.

In fact, the normal matrices are precisely those that are diagonalizable by a unitary matrix, and unitary matrices are perfectly conditioned in the 2-norm.

One minor problem with the diagonalization approach is that if A is real with some complex eigenvalues and f A is real, then rounding errors may cause the computed f A to have a tiny nonzero imaginary part. This problem can be overcome by discarding the computed imaginary part, but of course numerical instability may potentially produce a large spurious imaginary part.

For symmetric matrices the eigensystem is real and these considerations do not apply. This approach falls into the class of block diagonalization methods considered in Section 4. Another approach is described in 84 Techniques for General Functions the next section.

In summary, for a normal matrix, computing f A via a unitary diagonalization is the method of choice if the diagonalization can be computed. In particular, this method is recommended if A is symmetric or Hermitian. Schur Decomposition and Triangular Matrices We saw in the last section the importance of using well conditioned similarity transformations in order to maintain numerical stability.

Taking this view to the extreme, we can restrict to unitary transformations. In general, the closest one can go to diagonal form via unitary similarities is the Schur triangular form. The eigenvalues of A appear on the diagonal of T. The Schur decomposition can be computed with perfect backward stability by the QR algorithm, and hence it is a standard tool in numerical linear algebra. We therefore now turn our attention to functions of triangular matrices.

In fact, explicit formulae are available for all the elements of f T. Each of these products is zero. For general f , Theorem 3. An interesting property revealed by the theorem is that the 1,2 block of f A depends only linearly on the 1,2 block of A. To compute f T via Theorem 4. He notes that f T commutes with T see Theorem 1. Hence this recurrence enables F to be computed either a superdiagonal at a time, starting with the diagonal, or a column at a time, from the second column to the last, moving up each column.

In this situation 4. We will assume that T is also triangular, though this is not necessary to derive the recurrence. The Sylvester equation 4. Reordering can be achieved by standard techniques and we return to this topic in Section 9. When f A is real, this enables it to be computed entirely in real arithmetic.

We turn now to numerical considerations. In Listing 4. The function is in principle applicable to any matrix with distinct eigenvalues. It is very similar to the function funm in versions 6. Function funm simple often works well. With f the exponential, Table 4. The condition number of f A see Chapter 3 is about 2 in each case, so we would expect to be able to compute f A accurately. For A itself, funm simple yields an error of order 1, which is expected since it is unable to compute any of the superdiagonal elements of f A.

For the third matrix, in which the perturbation is triangular and the eigenvalues are 88 Techniques for General Functions Listing 4.

The conclusion of this experiment is that it is not just repeated eigenvalues that are bad for the scalar Parlett recurrence. Close eigenvalues can lead to a severe loss of accuracy, and even when the eigenvalues are far apart the recurrence can produce more inaccurate answers than expected.