Later on in the course, when we deal with the problem of computing the price of an European Option efficiently using the method of Dynamic Programming, we will have to deal with the problem of taking an n × n matrix A and raising it to a large number. That is, we will have to compute Ak for large values of k.
You could compute it by multiplying A over with it self k times. That is,
(with NEWMAT) something along the lines of
1: C = A;
2: for int i = 1; i < k; i++ do
3: C = A*C;
4: end for
Assuming you are doing straightforward matrix multiplication (nothing clever like Strassen’s method, etc.) then the above procedure will take O(n3k) steps, as each matrix multiplication is O(n3) and there are k-many of them.
A candidate algorithm that is more efficient than the one shown above goes by the name of Repeated Squaring, and it described below.
Repeated Squaring Algorithm Suppose you want to compute A11, you write the exponent 11 in binary – which is h1 0 1 1i2. That is, 11 = 1×23 +0×22 +1×21 +1×20. You then compute all dlog2 11e-many powers of A. That is, you compute A,A2,A4,A8, then add the appropriate components as per the binary expansion of the exponent. In our case, we add A8 ×A2 ×A to get the final product. It turns out that these constituent powers of A (i.e. 8, 2 and 1) can be computed in a recursive manner quite efficiently. Here is something from the following link that computes nk for integers n and k.
1: int power( int n, int k )
2: if k == 0 then
3: return 1
4: end if
5: if k is odd then
6: return (n * power (
7: else
8: return (power(
9: end if
A similar approach can be used to compute Ak for a square matrix A. I leave the nitty-gritty details for you to figure out. Keep in mind that with NEWMAT on your side, the matrix operations are structurally the same as that for scalars.
Complexity of Repeated Squaring vs. Brute Force Multiplication If you have to compute Ak, and assuming each multiplication is O(b3) (nothing clever here, straightforward matrix multiplication), and since we have log2 kmany of these to do (or, since log we can replace log2 k with
O(logk)), this entire operation takes O(b3 logk). This is in contrast to the O(b3k) procedure for multiplying A with itself k times. Resulting in a total complexity of O(b3 lnk). If we used Strassen’s method for matrix multiplication we would have a procedure that is O(b2.8[1] lnk).
The Programming Assignment 1. (Using NEWMAT) I want you to write a recursion routine
Matrix repeated squaring(Matrix A, int exponent, int no rows)
which takes a (no rows × no rows) matrix A and computes Aexponent using the Repeated Squaring Algorithm.
2. Your code should be able to take as input the size and exponent as input on the command line. That is, if we want to compute Ak, where A is an (n×n) square matrix, I want to be able to read n and k on the commandline. It should fill the entries of the matrix A with random entries in the interval (−5,5)1.
3. The output should indicate: (1) The number of rows/columns in A (that is read from the command line), (2) The exponent k (that is read from the command line), (3) The result and the amount of time it took to compute Ak using repeated squaring, and (4) The result and the amount of time it took to compute Ak using brute force multiplication. A sample output is shown in figure 1.
4. I want you to provide a plot of the computation time (in seconds) for the two methods as a function of the size of the matrix. That is, I am looking for something along the lines of figure 2. For this, you will have to place timer objects before and after appropriate portions of your code and do the needful – as the following lines of code illustrate.
time before = clock();
B = repeated squaring(A, exponent, dimension); time after = clock(); diff = ((float) time after - (float) time before);
cout << ”It took ” << diff/CLOCKS PER SEC << ” seconds to complete” << endl;
In terms of the value of the exponent, tell me the regions where one algorithm performs better than the other, as far as computation time is concerned.
Figure 1: A sample output for different exponents and matrix-dimensions.
Figure 2: A comparison of the computation-time (obtained experimentally) for brute-force exponentiation and the method of repeated squares of a random 5 × 5 matrix as a function of the exponent.