In [1]:
A = randn(3,3)
Out[1]:
and a vector of ones
In [2]:
x = ones(3)
Out[2]:
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
In [3]:
b = A*x
Out[3]:
In [4]:
Asym = A + A'
Out[4]:
In [5]:
Apd = A'A
Out[5]:
In [6]:
A\b
Out[6]:
In [7]:
A[:,1:2]\b
Out[7]:
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]:
In [9]:
A[1:2,:]\b[1:2]
Out[9]:
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.
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]:
The different parts of the factorization can be extracted by special indexing
In [11]:
Alu[:P]
Out[11]:
In [12]:
Alu[:L]
Out[12]:
In [13]:
Alu[:U]
Out[13]:
We can therefore compute the solution of $Ax=b$ from the factorization
In [14]:
Alu[:U]\(Alu[:L]\(Alu[:P]'b))
Out[14]:
However, more importantly the LU type allows dispatch and we can solve the system by
In [15]:
Alu\b
Out[15]:
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]:
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]:
\ has a method for the QR and the least squares problem is therefore solved with
In [18]:
Aqr\b
Out[18]:
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]:
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]:
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]:
The matrix has a compact internal representation and indexing is therefore not defined
In [22]:
Aqr[:Q][1]
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]:
and
In [24]:
Aqr[:Q]*ones(3)
Out[24]:
works, but not
In [26]:
Aqr[:Q]*ones(4)
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]:
The values and the vectors can be extracted from the Eigen type by special indexing
In [28]:
AsymEig[:values]
Out[28]:
In [29]:
AsymEig[:vectors]
Out[29]:
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]:
Julia also has an eig function which returns a tuple with the values and the vectors
In [31]:
eig(Asym)
Out[31]:
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]:
and again \ has a method for the type enabling least squares by SVD
In [33]:
Asvd\b
Out[33]:
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.
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]:
In [35]:
Triangular(tril(A))
Out[35]:
In [36]:
Symmetric(Asym)
Out[36]:
In [37]:
SymTridiagonal(diag(Asym),diag(Asym,1))
Out[37]:
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.
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))
In [38]:
@time eigvals(Asym1);
In [38]:
@time eigvals(Asym2);
In [38]:
@time eigvals(Symmetric(Asym2));
In [40]:
n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)
Out[40]:
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]:
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.
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]:
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]:
In [44]:
nHilbert = 8
H = Rational{BigInt}[1//(i+j-1) for i = 1:nHilbert,j = 1:nHilbert]
Out[44]:
In [45]:
inv(H)
Out[45]:
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]:
In [47]:
Hbf = convert(Matrix{BigFloat}, H);
float64(norm(inv(Hbf)-inv(H),Inf))
Out[47]:
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]:
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)
We can now create a square matrix and a vector over $GF(2)$
In [142]:
Agf2 = convert(Matrix{GF2},randbool(4,4))
Out[142]:
In [143]:
bgf2 = convert(Vector{GF2}, randbool(4))
Out[143]:
and because the LU factorization in Julia works whenever the
In [144]:
xgf2 = Agf2\bgf2
Out[144]:
In [145]:
Agf2*xgf2 - bgf2
Out[145]:
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.
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]:
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]:
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]:
The type does not print nicely and it is slightly more informative to get the property names
In [153]:
names(AherFact)
Out[153]:
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]:
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.
In [154]:
n = 500;
A = randn(n,n);
Apd = A'A;
@time eigvals(Symmetric(Apd));
In [154]:
A = complex(randn(n,n),randn(n,n));
Apd = A'A;
@time eigvals(Hermitian(Apd));
In [154]:
A = Quaternion{Float64}[Quaternion(randn(4)...) for i = 1:n,j=1:n];
Apd = A'A;
@time eigvals(Hermitian(Apd,:L));
In [ ]: