194 6 Linear Algebra
(bˆ−Aˆˆxk)−αkAˆpˆk,which gives
rˆk+ 1 =rˆk−Aˆpˆk,
If we consider finding the minimum of a functionfusing Newton’s method, that implies a
search for a zero of the gradient of a function. Near a pointxiwe have to second order
f(ˆx) =f(xˆi)+ (xˆ−xˆi)∇f(xˆi)1
2 (
xˆ−ˆxi)Aˆ(xˆ−ˆxi)giving
∇f(ˆx) =∇f(xˆi)+Aˆ(xˆ−xˆi).
In Newton’s method we set∇f= 0 and we can thus compute the next iteration point
xˆ−ˆxi=Aˆ−^1 ∇f(xˆi).Subtracting this equation from that ofxˆi+ 1 we have
xˆi+ 1 −ˆxi=Aˆ−^1 (∇f(xˆi+ 1 )−∇f(xˆi)).6.7 A vector and matrix class
We end this chapter by presenting a class which allows to manipulate one- and two-
dimensional arrays. However, before we proceed, we would like to come with some general
recommendations. Although it is useful to write your own classes, like the one included here,
in general these classes may not be very efficient from a computational point of view. There
are several libraries which include many interesting arrayfeatures that allow us to write
more compact code. The latter has the advantage that the codeis lost likely easier to debug
in case of errors (obviously assuming that the library is functioning correctly). Furthermore,
if the proper functionalities are included, the final code may closely resemble the mathemat-
ical operations we wish to perform, increasing considerably the readability of our program.
And finally, the code is in almost all casesmuch faster than the one we wrote!
In particular, we would like to recommend the C++ linear algebra library Armadillo, see
http://arma.sourceforgenet. For those of you who are familiar with compiled programs
like Matlab, the syntax is deliberately similar. Integer, floating point and complex numbers
are supported, as well as a subset of trigonometric and statistics functions. Various matrix
decompositions are provided through optional integrationwith LAPACK, or one of its high
performance drop-in replacements (such as the multi-threaded MKL or ACML libraries). The
selected examples included here show some examples on how todeclare arrays and rearrange
arrays or perform mathematical operations on say vectors ormatrices. The first example
here defines two random matrices of dimensionality 10 × 10 and performs a matrix-matrix
multiplication using thedgemmfunction of the library BLAS.
Simple matrix-matrix multiplication of two random matrices
#include
#include
using namespacestd;
using namespacearma;
intmain(intargc,char**argv)
{