Linear algebra in Julia

Outline

  • Basic linear algebra
  • Factorizations
  • Special matrices
  • Generic linear algebra

Basic linear algebra

Syntax very similar to MATLAB by there are some important differences. Define a matrix of random normal variates


In [1]:
A = randn(3,3)


Out[1]:
3x3 Array{Float64,2}:
 1.22112  -1.04948   -0.544682
 2.00471  -0.121917  -1.62091 
 1.64939  -0.681989   0.362995

and a vector of ones


In [2]:
x = ones(3)


Out[2]:
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

Notice that $A$ has type Array{Float64,2} but $x$ has type Array{Float64,1}. Julia defines the aliases Vector{Type}=Array{Type,1} and Matrix{Type}=Array{Type,2} in Julia. The destinction between between vectors and matrices typically causes some confusion when coming from a MATLAB background.

Many of the basic operations are the same as in MATLAB

Multiplication


In [3]:
b = A*x


Out[3]:
3-element Array{Float64,1}:
 -0.373049
  0.261885
  1.3304  

Transposition

As in MATLAB A' is the conjugate transpose whereas A.' is just the transpose


In [4]:
Asym = A + A'


Out[4]:
3x3 Array{Float64,2}:
 2.44224    0.955227   1.10471 
 0.955227  -0.243835  -2.3029  
 1.10471   -2.3029     0.725991

Transposed multiplication

Julia allows us to write this without *


In [5]:
Apd = A'A


Out[5]:
3x3 Array{Float64,2}:
  8.23049  -2.65082   -3.31585 
 -2.65082   1.58139    0.521694
 -3.31585   0.521694   3.05579 

Solving linear systems

The problem $Ax=b$ for square $A$ is solved by the \ function.


In [6]:
A\b


Out[6]:
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

Least squares

When A is tall the \ function calculates the least squares solution


In [7]:
A[:,1:2]\b


Out[7]:
2-element Array{Float64,1}:
 0.355336
 0.249271

The \ function also works for rank deficient least squares problems. In this case, the least squares solution is not unique and Julia returns the solution with smallest norm


In [8]:
[A[:,1] A[:,1]]\b


Out[8]:
2-element Array{Float64,1}:
 0.137526
 0.137526

Underdetermined systems

Minimum norm solution is returned


In [9]:
A[1:2,:]\b[1:2]


Out[9]:
3-element Array{Float64,1}:
 -0.0050698
  0.454406 
 -0.202015 

Factorization

The \ function hides how the problem is actually solved. Depending on the dimensions of A, different methods are chosen to solve the problem. An intermediate step in the solution is to calculate a factorization of the matrix A. Basically, a factorization of A is a way of expressing A as a product of triangular, unitary and permutation matrices. Julia defines a Factorization abstract type and several composite subtypes for actually storing factorizations. A Factorization object should be thought of as a representation of the matrix A.

LU

When A is square, problem is solved by factorizing the matrix A=PLU where P is a permutation matrix, L is lower triangular unit diagonal and U is upper triangular. Julia allows computing the LU factorization and defines a composite factorization type for storing it.


In [10]:
Alu = lufact(A)


Out[10]:
LU{Float64}(3x3 Array{Float64,2}:
 2.00471   -0.121917  -1.62091 
 0.609124  -0.975222   0.442652
 0.822758   0.59646    1.43259 ,[2,2,3],0)

The different parts of the factorization can be extracted by special indexing


In [11]:
Alu[:P]


Out[11]:
3x3 Array{Float64,2}:
 0.0  1.0  0.0
 1.0  0.0  0.0
 0.0  0.0  1.0

In [12]:
Alu[:L]


Out[12]:
3x3 Array{Float64,2}:
 1.0       0.0      0.0
 0.609124  1.0      0.0
 0.822758  0.59646  1.0

In [13]:
Alu[:U]


Out[13]:
3x3 Array{Float64,2}:
 2.00471  -0.121917  -1.62091 
 0.0      -0.975222   0.442652
 0.0       0.0        1.43259 

We can therefore compute the solution of $Ax=b$ from the factorization


In [14]:
Alu[:U]\(Alu[:L]\(Alu[:P]'b))


Out[14]:
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

However, more importantly the LU type allows dispatch and we can solve the system by


In [15]:
Alu\b


Out[15]:
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

This could be useful if the same left-hand-side is used for several right-hand-sides. The factorization can also be used for calculating the determinant because $\det(A)=\det(PLU)=\det(P)\det(U)=\pm \prod u_{ii}$ because $U$ is triangular and the sign is determined from $\det(P)$.


In [16]:
det(Alu)


Out[16]:
2.800760493907129

QR

When A is tall, Julia computes the least squares solution $\hat{x}$ that minimizes $\|Ax-b\|_2$. This can be done by factorizing $A=QR$ where $Q$ is unitary/orthogonal and $R=\left(\begin{smallmatrix}R_0\\0\end{smallmatrix}\right)$ and $R_0$ is upper triangular. With the QR factorization the minimum norm can be expressed \begin{equation*} \|Ax-b\|=\|QRx-b\|=\|Q(Rx-Q'b)\|=\|Rx-Q'b\|=\left\|\begin{pmatrix}R_0x-Q_0'b\\Q_1'b\end{pmatrix}\right\|=\|R_0x-Q_0'b\|+\|Q_1'b\| \end{equation*} and the problem therefore reduces to solving the square problem $R_0x=Q_0'b$ for $x$.

We can QR factorize the submatrix of the two first columns of $A$ by


In [17]:
Aqr = qrfact(A[:,1:2])


Out[17]:
QRCompactWY{Float64}(3x2 Array{Float64,2}:
 -2.86888    0.923991 
  0.490149  -0.853013 
  0.403274   0.0670412,2x2 Array{Float64,2}:
 1.42564       -1.46804
 5.45353e-312   1.99105)

\ has a method for the QR and the least squares problem is therefore solved with


In [18]:
Aqr\b


Out[18]:
2-element Array{Float64,1}:
 0.355336
 0.249271

It should be noted that this is not the way A[:,1:2]\b is solved. In order to handle rank deficient problems Julia uses a QR factorization with pivoting. Pivoting is enabled with a keyword


In [19]:
Aqrp = qrfact([A[:,1] A[:,1]],pivot=true)


Out[19]:
QRPivoted{Float64}(3x2 Array{Float64,2}:
 -2.86888   -2.86888
  0.490149   0.0    
  0.403274   0.0    ,[1.42564,0.0],[1,2])

Notice that the type is different now. \ also has a method for QRPivoted and the rank deficient problem is therefore computed


In [20]:
Aqrp\b


Out[20]:
2-element Array{Float64,1}:
 0.137526
 0.137526

Another feature of the QR factorizations is the Q types for storing the unitary matrices $Q$. They can be extracted from the different from QR types by indexing


In [21]:
Aqr[:Q]


Out[21]:
3x3 QRCompactWYQ{Float64}:
 -0.425642   0.769267
 -0.698778  -0.613996
 -0.574925   0.176743

The matrix has a compact internal representation and indexing is therefore not defined


In [22]:
Aqr[:Q][1]


indexing not defined for QRCompactWYQ{Float64}
while loading In[22], in expression starting on line 1
 in getindex at abstractarray.jl:357

Even though the matrix is printed as a $3\times 2$ matrix is a in practice representing the square version as well. Hence both


In [23]:
Aqr[:Q]*ones(2)


Out[23]:
3-element Array{Float64,1}:
  0.343624
 -1.31277 
 -0.398182

and


In [24]:
Aqr[:Q]*ones(3)


Out[24]:
3-element Array{Float64,1}:
 -0.132881
 -1.67981 
  0.400707

works, but not


In [26]:
Aqr[:Q]*ones(4)


DimensionMismatch("")
while loading In[26], in expression starting on line 1
 in * at linalg/factorization.jl:387

Eigendecompositions and the SVD(s)

The results from eigendecompositions and singular values decompositions are also stored in Factorization types. This also includes Hessenberg and Schur factorizations.

The eigendecomposition can be computed


In [27]:
AsymEig = eigfact(Asym)


Out[27]:
Eigen{Float64,Float64}([-2.52923,2.44803,3.00559],3x3 Array{Float64,2}:
  0.278461   0.490567   0.825714
 -0.736987   0.660426  -0.143829
 -0.615881  -0.56849    0.545445)

The values and the vectors can be extracted from the Eigen type by special indexing


In [28]:
AsymEig[:values]


Out[28]:
3-element Array{Float64,1}:
 -2.52923
  2.44803
  3.00559

In [29]:
AsymEig[:vectors]


Out[29]:
3x3 Array{Float64,2}:
  0.278461   0.490567   0.825714
 -0.736987   0.660426  -0.143829
 -0.615881  -0.56849    0.545445

Once again, when the factorization is stored in a type, we can dispatch on it and write specialized methods that exploit the properties of the factorization, e.g. that $A^{-1}=(V\Lambda V^{-1})^{-1}=V\Lambda^{-1}V^{-1}$.


In [30]:
inv(AsymEig)*Asym


Out[30]:
3x3 Array{Float64,2}:
  1.0          -5.27356e-16  -5.13478e-16
 -9.4369e-16    1.0           1.33227e-15
 -3.33067e-16   4.44089e-16   1.0        

Julia also has an eig function which returns a tuple with the values and the vectors


In [31]:
eig(Asym)


Out[31]:
([-2.52923,2.44803,3.00559],
3x3 Array{Float64,2}:
  0.278461   0.490567   0.825714
 -0.736987   0.660426  -0.143829
 -0.615881  -0.56849    0.545445)

This is only provided for MATLAB compatibility and I do not recomment this version.

The svdfact function computes the singular value decomposition


In [32]:
Asvd = svdfact(A[:,1:2])


Out[32]:
SVD{Float64,Float64}(3x2 Array{Float64,2}:
 -0.495403   0.726305
 -0.638588  -0.676377
 -0.588881   0.122457,[3.02621,0.808667],2x2 Array{Float64,2}:
 -0.943896   0.330243
 -0.330243  -0.943896)

and again \ has a method for the type enabling least squares by SVD


In [33]:
Asvd\b


Out[33]:
2-element Array{Float64,1}:
 0.355336
 0.249271

In contrast to MATLAB, Julia does not allow dispatch on the number of output arguments and therefore there are special functions for providing values only: eigvals and svdvals.

Special matrices

The structure of matrices is very important in linear algebra. This structure can be made explicit in Julia through composite types. Examples are Diagonal, Triangular, Symmetric, Hermitian, Tridiagonal and SymTridiagonal. Specialized methods are written for the special matrix types to take advantage of their structure. Below some examples are shown


In [34]:
Diagonal(diag(A))


Out[34]:
3x3 Diagonal{Float64}:
 1.22112   0.0       0.0     
 0.0      -0.121917  0.0     
 0.0       0.0       0.362995

In [35]:
Triangular(tril(A))


Out[35]:
3x3 Triangular{Float64}:
 1.22112   0.0       0.0     
 2.00471  -0.121917  0.0     
 1.64939  -0.681989  0.362995

In [36]:
Symmetric(Asym)


Out[36]:
3x3 Symmetric{Float64}:
 2.44224    0.955227   1.10471 
 0.955227  -0.243835  -2.3029  
 1.10471   -2.3029     0.725991

In [37]:
SymTridiagonal(diag(Asym),diag(Asym,1))


Out[37]:
3x3 SymTridiagonal{Float64}:
 2.44224    0.955227   0.0     
 0.955227  -0.243835  -2.3029  
 0.0       -2.3029     0.725991

When it is known that a matrix is e.g. triangular or symmetric Julia might be able to solve a problem faster by converting the matrix to a special matrix. For some of the procedures, Julia checks if the input matrix is triangular or symmetric and converts the matrix if such a structure is detected. It should noted that Symmetric, Hermitian and Triangular do not copy the input matrix.

Symmetric eigenproblem

Whether or not Julia is able to detech if a matrix i symmetric/Hermitian can have a big influence on how fast an eigenvalue problem is solved. Sometimes it is known that a matric is symmetric or Hermitian but due to floating point errors this is not detected by the eigvals function. In following example Asym1 and Asym2 are almost identical, but the unless Julia is told that Asym2 is symmetric, the elapsed time for the computation is very different.


In [38]:
n = 1000;
A = randn(n,n);
Asym1 = A + A';
Asym2 = copy(Asym1);Asym2[1,2] += 5eps();
println("Is Asym1 symmetric? ", issym(Asym1))
println("Is Asym2 symmetric? ", issym(Asym2))


Is Asym1 symmetric? true
Is Asym2 symmetric? false

In [38]:
@time eigvals(Asym1);


elapsed time: 0.470104509 seconds (8921768 bytes allocated)

In [38]:
@time eigvals(Asym2);


elapsed time: 2.61170262 seconds (8574096 bytes allocated)

In [38]:
@time eigvals(Symmetric(Asym2));


elapsed time: 0.428100059 seconds (8435512 bytes allocated)

A big problem

Using the tridiagonal matrices makes it possible to work with potentially very large problems. The following problem would not be possible to solve on a laptop if the matrix had to be stored as a Matrix type.


In [40]:
n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)


elapsed time: 1.003869886 seconds (216002736 bytes allocated)
Out[40]:
6.501819823351534

Generic linear algebra

The usual way of adding support for numerical linear algebra is by wrapping BLAS and LAPACK subroutines. For matrices with elements of Float32, Float64, Complex{Float32} or Complex{Float64} this is also what Julia do. For a long time Julia has also had support for multiplicaton of general element types. Hence, when multiplying integer matrices, the output is also an integer matrix


In [41]:
rand(1:10,3,3)*rand(1:10,3,3)


Out[41]:
3x3 Array{Int64,2}:
  96   80  117
 125   62  122
 139  103  155

Recently, more generic linear algebra methods has been added and Julia now supports generic LU and QR factorizations. Generic eigenvalue methods are being tested and can hopefully be added to Julia soon.

In general, the LU factorization can be computed whenever the matrix element type is closed under the operations +, -, * and \. Of couse the matrix also has to have full rank. The generic LU method in Julia applies pivoting and therefore the element type also has to support < and abs. In consequence it is possible to solve systems of equations of e.g. rational numbers which the following examples show.

Example 1: Rational linear system of equations

Julia has rational numbers build in. The following example shows how a linear system of equations can be solved without promoting to floating point element types. Overflow can easisly become a problem when working with rational numbers so we use BigInts.


In [42]:
Ar = convert(Matrix{Rational{BigInt}}, rand(1:10,3,3))/10
x = ones(Int,3)
b = Ar*x
Ar\b


Out[42]:
3-element Array{Rational{BigInt},1}:
 1//1
 1//1
 1//1

Example 2: Rational matrix from eigenstructure

The next example shows how we rational matrix arithmetic can be used for calculating a matrix given rational eigenvalues and -vectors. I have found this convenient when giving examples of linear dynamic systems.


In [43]:
λ1,λ2,λ3=1//1,1//2,1//4
v1,v2,v3=[1,0,0],[1,1,0],[1,1,1]
V,Λ=[v1 v2 v3], Diagonal([λ1,λ2,λ3])
A = V*Λ/V


Out[43]:
3x3 Array{Rational{Int64},2}:
 1//1  -1//2  -1//4
 0//1   1//2  -1//4
 0//1   0//1   1//4

Example 3: The Hilbert matrix.

The Hilbert matrix $H$ defined by $(h_{ij})=\frac{1}{1-i-j}$ is very ill conditioned. Even for small dimensions, the inverse of the Hilbert matrix is very imprecise in double precision. By defining the matrix with rational elements, the inverse can be calculated exact


In [44]:
nHilbert = 8
H = Rational{BigInt}[1//(i+j-1) for i = 1:nHilbert,j = 1:nHilbert]


Out[44]:
8x8 Array{Rational{BigInt},2}:
 1//1  1//2  1//3   1//4   1//5   1//6   1//7   1//8 
 1//2  1//3  1//4   1//5   1//6   1//7   1//8   1//9 
 1//3  1//4  1//5   1//6   1//7   1//8   1//9   1//10
 1//4  1//5  1//6   1//7   1//8   1//9   1//10  1//11
 1//5  1//6  1//7   1//8   1//9   1//10  1//11  1//12
 1//6  1//7  1//8   1//9   1//10  1//11  1//12  1//13
 1//7  1//8  1//9   1//10  1//11  1//12  1//13  1//14
 1//8  1//9  1//10  1//11  1//12  1//13  1//14  1//15

In [45]:
inv(H)


Out[45]:
8x8 Array{Rational{BigInt},2}:
      64//1      -2016//1       20160//1  …       192192//1      -51480//1
   -2016//1      84672//1     -952560//1       -10594584//1     2882880//1
   20160//1    -952560//1    11430720//1       141261120//1   -38918880//1
  -92400//1    4656960//1   -58212000//1      -776936160//1   216216000//1
  221760//1  -11642400//1   149688000//1      2118916800//1  -594594000//1
 -288288//1   15567552//1  -204324120//1  …  -3030051024//1   856215360//1
  192192//1  -10594584//1   141261120//1      2175421248//1  -618377760//1
  -51480//1    2882880//1   -38918880//1      -618377760//1   176679360//1

Even for $n=8$ the difference between the exact and the floating point result is big.


In [46]:
norm(float64(float(inv(H))) - inv(float64(float(H))),Inf)


Out[46]:
124.85434710979462

Another approach could be to use higher precision floating point numbers.

Example 4: A Hilbert matrix with BigFloat elements

Higher precision floating point numbers are build into Julia and it is therefore easy to convert the Hilbert matrix


In [47]:
Hbf = convert(Matrix{BigFloat}, H);
float64(norm(inv(Hbf)-inv(H),Inf))


Out[47]:
6.304476724749898e-59

Which is fairly precise. However, the matrix need not grow much before 256 bit floats are not sufficient


In [48]:
nHilbert = 30
H = Rational{BigInt}[1//(i+j-1) for i = 1:nHilbert,j = 1:nHilbert];
Hbf = convert(Matrix{BigFloat}, H);
float64(norm(inv(Hbf)-inv(H),Inf))


Out[48]:
2.1751007413310495e8

The next example shows how the LU factorization can be used for a user defined number type.

Example 5: GF(2)

Julia makes it easy to define new number types. The block below defines a type representing the Galois field with two elements GF(2).


In [48]:
import Base: &,|,<,+,-,*,/,abs,and_int,box,convert,div,inv,or_int,rem,print,promote_rule,show,showcompact,udiv_int,ult_int,unbox,urem_int,xor_int
bitstype 8 GF2 <: Unsigned
convert(::Type{GF2}, x::GF2) = x
convert(::Type{GF2}, x::Uint8) = x == 0 || x == 1 ? box(GF2,unbox(Uint8,x)) : throw(ArgumentError("input must be zero or one"))
convert(::Type{GF2}, x::Integer) = convert(GF2,uint8(x))
convert(::Type{Int8}, x::GF2) = box(Int8,unbox(GF2,x))
convert(::Type{Uint8}, x::GF2) = box(Uint8,unbox(GF2,x))
convert(::Type{Uint64}, x::GF2) = uint64(uint8(x))
gf2(x::Integer) = convert(GF2,x)
promote_rule(::Type{GF2},::Type{Uint64}) = Uint64
print(io::IO, x::GF2) = print(io, int8(x))
show(io::IO, x::GF2) = show(io, int8(x))
showcompact(io::IO, x::GF2) = showcompact(io, int8(x))
+(x::GF2, y::GF2) = box(GF2,xor_int(unbox(GF2,x),unbox(GF2,y)))
-(x::GF2, y::GF2) = box(GF2,xor_int(unbox(GF2,x),unbox(GF2,y)))
*(x::GF2, y::GF2) = box(GF2,and_int(unbox(GF2,x),unbox(GF2,y)))
inv(x::GF2) = x==one(GF2) ? x : error("integer division error")
/(x::GF2,y::GF2) = x*inv(y)
-(x::GF2) = x
abs(x::GF2) = x
div(x::GF2, y::GF2) = box(GF2,udiv_int(unbox(GF2,x),unbox(GF2,y)))
rem(x::GF2, y::GF2) = box(GF2,urem_int(unbox(GF2,x),unbox(GF2,y)))
<(x::GF2, y::GF2)   = ult_int(unbox(GF2,x),unbox(GF2,y))
(&)(x::GF2, y::GF2) = box(GF2,and_int(unbox(GF2,x),unbox(GF2,y)))
(|)(x::GF2, y::GF2) = box(GF2,or_int(unbox(GF2,x),unbox(GF2,y)));

The definitions < and abs do not really make sense for $GF(2)$, but they make it possible to do pivoting on the non-zero element. Below, the arithmetic rules are shown.


In [138]:
x = gf2(1)
y = gf2(0)
println("x+y: ", x+y)
println("x+x: ", x+x)
println("x-x: ", x-x)
println("x*y: ", x*y)
println("x*x: ", x*x)
println("x/x: ", x/x)
println("y/x: ", y/x)


x+y: 1
x+x: 0
x-x: 0
x*y: 0
x*x: 1
x/x: 1
y/x: 0

We can now create a square matrix and a vector over $GF(2)$


In [142]:
Agf2 = convert(Matrix{GF2},randbool(4,4))


Out[142]:
4x4 Array{GF2,2}:
 0  0  0  1
 1  1  1  0
 0  1  1  1
 1  0  1  1

In [143]:
bgf2 = convert(Vector{GF2}, randbool(4))


Out[143]:
4-element Array{GF2,1}:
 1
 0
 0
 1

and because the LU factorization in Julia works whenever the


In [144]:
xgf2 = Agf2\bgf2


Out[144]:
4-element Array{GF2,1}:
 1
 0
 1
 1

In [145]:
Agf2*xgf2 - bgf2


Out[145]:
4-element Array{GF2,1}:
 0
 0
 0
 0

It worked! The final example shows how the generic elementary reflectors work. The first step in solving a Hermitian eigenvalue problem is to tranform the Hermitian matrix to a real symmetric tridiagonal matrix by elementary reflectors of the form $H=I-v\tau v'$ satisfying $H'H=I$, i.e. they are unitary. Most often, textbooks present the special case of Householder reflectors $H=I-2vv'$ where $v'v=1$. Householder reflectors have the nice property that they are also Hermitian, but they fail to reduce Hermitian matrix to a real symmetric tridiagonal matrix. The following example shows how the underlying real symmtric tridiagonal eigenproblem can be extracted from a Hermitian eigenvalue problem of a matrix with a userdefined element type.

Example 6: Quaternions

In my github repo, I have a package that defines quaternions. Recall that quaternions generalize complex numbers and can be written $q=x_1+x_2i+x_3j+x_4k$ where the "imaginary" units $i$,$j$ and $k$ satisfies $i^2=j^2=k^2=ijk=-1$. The only difficulty in defining quaternions is the fact that they do not commute, e.g. $ij=-ijkk=k$ but $ji=-jiijk=jjk=-k$. However, the generic algorithms for LU factorization and elementary reflectors do not assume commutativity and therefore we can solve positive definite eigenproblems for quaternions. First we compute a Hermitian matrix of quaternions


In [146]:
using Quaternions
A = Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:2,j=1:2];
Aher = A + A'


Out[146]:
2x2 Array{Quaternion{Float64},2}:
                 -1.42901+0.0i+0.0j+0.0k  …  -1.20241+0.293595i-3.21715j+0.00250102k
 -1.20241-0.293595i+3.21715j-0.00250102k                     -2.27177+0.0i+0.0j+0.0k

We can now calculate the eigenvalues. In this case, we need explicitly to call the Hermtian constructor bacuse Julia Base only checks is the matrix is Hermitian when the element type is supported by LAPACK.


In [147]:
eigvals(Hermitian(Aher))


Out[147]:
2-element Array{Float64,1}:
 -5.32309
  1.62231

We can break down the steps of the calculation. First, the matrix is transformed to a symmtric tridiagonal matrix


In [152]:
AherFact = Quaternions.symtri!(copy(Hermitian(Aher)))


Out[152]:
SymmetricTridiagonalFactorization{Quaternion{Float64}}(2x2 Array{Quaternion{Float64},2}:
 -1.42901+0.0i+0.0j+0.0k  -1.20241+0.293595i-3.21715j+0.00250102k
  3.44704+0.0i+0.0j+0.0k                  -2.27177+0.0i+0.0j+0.0k,Quaternion{Float64}[1.34882+0.085173i-0.933309j+0.000725555k],2x2 SymTridiagonal{Float64}:
 -1.42901   3.44704
  3.44704  -2.27177)

The type does not print nicely and it is slightly more informative to get the property names


In [153]:
names(AherFact)


Out[153]:
3-element Array{Symbol,1}:
 :factors      
 :scalarfactors
 :diagonals    

The property factors is a matrix with quaternion elements whose diagonal and subdiagonal contain the real symmetric tridiagonal matrix, i.e. the $i$,$j$ and $k$ terms are all zero. The elements below the subdiagonal contains the elements of the vectors $v$ of the elementary reflectors. The property scalafactors contains the multiplyers $\tau$. Finally, the property diagonals contains a real SymTridiagonal matrix.

When we do not assume commutativity, an elementary reflector is written $H=I-v\tau v^*$. As in LAPACK we use the convention that the first element of $v$ is one. Hence, we can write check the factorization $A^*A=QSQ^*$ where $Q$ is unitary and $S$ is a tridiagonal matrix.


In [154]:
Q = [1 0;0 1-AherFact.scalarfactors[1]];
S = full(AherFact.diagonals);
norm(Aher-Q*S*Q',1)


Out[154]:
1.7766851532539086e-15

All arithmetic operations are more expensive because a quaternion contains four real elements. Furthermore, LAPACK is fast bacause it uses blocked algorithms which we do not in our generic code. Anyway, the reduction to a real symmetric tridiagonal matrix can be done fairly fast. Below we present some results for computational time and memory usage.

Real


In [154]:
n = 500;
A = randn(n,n);
Apd = A'A;
@time eigvals(Symmetric(Apd));


elapsed time: 0.082535236 seconds (2189472 bytes allocated)

Complex


In [154]:
A = complex(randn(n,n),randn(n,n));
Apd = A'A;
@time eigvals(Hermitian(Apd));


elapsed time: 0.188167424 seconds (4417584 bytes allocated)

Quaternion


In [154]:
A = Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:n,j=1:n];
Apd = A'A;
@time eigvals(Hermitian(Apd,:L));


elapsed time: 0.884585021 seconds (8199712 bytes allocated)

In [ ]: