In [1]:
v1 = collect(1:6);
size(v1) # size always returns a tuple
Out[1]:
In [2]:
length(v1)
Out[2]:
In [3]:
m1 = reshape(collect(1:6),(3,2))
Out[3]:
Because size returns a tuple the number of rows and columns can be assigned in one statement. A common idiom is
In [4]:
m,n = size(m1)
Out[4]:
In [5]:
println("Number of rows = $m and number of columns = $n")
As in R, matrices and higher-dimensional arrays are stored in column-major ordering. Transpose is indicated by '. The * operator applied to arrays denotes matrix multiplication. Element-wise multiplication is indicated by .*. In general a prefix . to an operator makes it behave element-wise. Constructors of arrays with specific content include ones(), zeros() and eye() (the identity matrix).
In [6]:
m1 * ones(2)
Out[6]:
In [7]:
m1'*m1
Out[7]:
In [8]:
eye(3)
Out[8]:
In [9]:
eye(3,4)
Out[9]:
In [10]:
ones(Int,(2,3))
Out[10]:
The equivalent of R's cbind and rbind functions are hcat and vcat.
In [11]:
m2 = hcat(ones(3),collect(-1.:1.))
Out[11]:
vcat is used for concatenating vectors.
In [12]:
vcat([1:3;],ones(2))
Out[12]:
In some circumstances you can omit the * in a multiplication.
In [13]:
m22 = m2'm2
Out[13]:
Unlike in R, the Julia expression X'X is evaluated using only one copy of X and will always produce a symmetric matrix. In R the function crossprod must be called to do this. (Be careful, in Julia there is an operation called a cross product but it is the cross product from Physics.)
In [14]:
?×
Out[14]:
The backslash operator, \, denotes the solution of a linear system of equations of the form $Ax=b$.
In [15]:
m22\(m2'ones(3))
Out[15]:
With a rectangular left-hand side \ denotes a least-squares solution
In [16]:
m2\ones(3)
Out[16]:
Note that this solution is a bit different from the previous solution because of numerical imprecision in the computation.
Those who haven't studied numerical linear algebra assume that a linear system of equations is solved as $A^{-1}b$. In practice it is almost never necessary to evaluate the inverse of a matrix. Instead the matrix is factored into a product of other matrices that have desirable properties such as being triangular.
Such factorizations include the singular value decomposition (SVD), the QR decomposition (with or without pivoting), the Cholesky factorization (with or without pivoting) for a positive-semidefinite symmetric matrix, the LU factorization and the Eigen decomposition. The names of the constructors of the decomposition usually end in fact.
In [17]:
m3 = hcat(ones(3),[1:3;])
Out[17]:
In [18]:
m33 = m3'm3
Out[18]:
In [19]:
cf33 = cholfact(m33)
Out[19]:
In [20]:
cf33[:L] # extract the lower Cholesky factor
Out[20]:
In [21]:
cf33[:U] # extract the upper triangular Cholesky factor, cf33[:L]'
Out[21]:
In [22]:
full(cf33) # the matrix (up to round-off) represented by the factorization
Out[22]:
In [23]:
det(cf33) # determinant of full(cf33), calculated from det(cf33[:U])
Out[23]:
In [24]:
logdet(cf33) # more practical for positive semidefinite
Out[24]:
Although one rarely needs the inverse of a general matrix, it is sometimes useful to get the inverse of the matrix represented by the Cholesky factorization. This is feasible because the factor is triangular and determining the inverse of a triangular matrix is easier than doing so for a general matrix. Suppose we want the estimated covariance matrix of the coefficient estimates.
In [25]:
y = rand(3)
Out[25]:
In [26]:
β = cf33\(m3'y)
Out[26]:
In [27]:
ŷ = m3*β # note you have to construct ŷ as y\hat<tab> not \hat<tab>y
Out[27]:
In [28]:
resid = y - ŷ
Out[28]:
In [29]:
s² = sumabs2(resid)/(3-2)
Out[29]:
In [30]:
vcv = s²*inv(cf33)
Out[30]:
The pivoted Cholesky decomposition has better numerical properties than does the unpivoted version and is a rank-revealing decomposition.
In [31]:
cf33 = cholfact(m33, :U, Val{true})
Out[31]:
At present this factorization doesn't have an attractive show method.
In [32]:
cf33[:U]
Out[32]:
In [33]:
cf33[:p] # the pivot vector
Out[33]:
In [34]:
cf33[:P] # the pivot matrix
Out[34]:
In [35]:
rank(cf33)
Out[35]:
The QR decomposition is a direct (i.e. non-iterative) decomposition used, among other things, to solve least squares problems. It decomposes a rectangular X into an orthogonal matrix Q and an upper triangular matrix R. The matrix Q, which could be very large if the number of rows of X is large, is stored implicitly.
In [36]:
qr3 = qrfact(m3)
Out[36]:
In [37]:
qr3[:R]
Out[37]:
In [38]:
qr3[:Q] # the implicit representation
Out[38]:
The QR factorization can be used directly for least squares solutions.
In [39]:
β = qr3\y
Out[39]:
What are called the effects in an R lm object are Q'y
In [40]:
qr3[:Q]'y
Out[40]:
There is also a pivoted QR decomposition which is rank-revealing.
In [41]:
qr3 = qrfact(m3, Val{true})
Out[41]:
The singular value decomposition represents X as U*S*V' where U and V are orthogonal and S is diagonal. It is related to the eigen decomposition in that the eigen decomposition of X'X is V*S*S*V'. The eigenvalues of X'X are the squares of the singular values of X.
In [42]:
sv3 = svdfact(m3)
Out[42]:
In [43]:
sv3[:U]
Out[43]:
In [44]:
sv3[:V]
Out[44]:
In [45]:
sv3[:S] # vector of singular values
Out[45]:
In [46]:
sv3[:Vt]
Out[46]:
This is the thin SVD. If you want the full matrix U use
In [47]:
sv3 = svdfact(m3; thin=false)
Out[47]:
In [48]:
sv3[:U]
Out[48]:
The eigfact function has a special method for real symmetric matrices
In [49]:
ef3 = eigfact(m33)
Out[49]:
In [50]:
ef3[:values]
Out[50]:
In [51]:
abs2(sv3[:S]) # same as the eigenvalues but in the opposite order
Out[51]:
In [52]:
ef3[:vectors]
Out[52]:
In [53]:
sv3[:V] # same as the eigenvectors of m33 with permuted rows and columns
Out[53]:
For general matrices we must allow for complex eigenvalues and eigenvectors. If all the eigenvalues and eigenvectors are real they are returned as such
In [54]:
eigfact(reshape([1:9;],(3,3)))
Out[54]:
In [55]:
eigfact([0 1;-1 0])
Out[55]:
The LU factorization is the most general direct factorization. It is related to Gaussian elimination.
In [56]:
lu4 = lufact(reshape([1:9;],(3,3)))
Out[56]:
In [57]:
lu4[:L] # unit lower triangular factor
Out[57]:
In [58]:
lu4[:U] # upper triangular factor. Note that it is singular.
Out[58]:
In [59]:
lu4[:p] # permutation vector
Out[59]:
In [60]:
lu4[:P] # permutation matrix
Out[60]:
As in R elements or subarrays are extracted with using indices. To indicate all values of a given index use :
In [61]:
m2[:,2]
Out[61]:
There are two alternatives to indexing, the use of sub or slice.
In [62]:
sub(m2,:,2)
Out[62]:
In [63]:
slice(m2, :, 2)
Out[63]:
In [64]:
sub(m2, 2, :)
Out[64]:
In [65]:
slice(m2, 2, :)
Out[65]:
The difference between sub and slice is that slice always drops singleton dimensions whereas sub only drops trailing singleton dimensions.
Both sub and slice return views into the original array. That is, they do not copy. At present using indices produces a copy but that will change.
In [66]:
pointer(m2)
Out[66]:
In [67]:
pointer(m2[:,1])
Out[67]:
In [68]:
pointer(sub(m2, :, 1))
Out[68]:
In [69]:
pointer(slice(m2, 1, :))
Out[69]:
To convert a matrix or higher-order array to a vector use vec
In [70]:
vec(m2)
Out[70]:
The result of an indexing operation is a copy of the elements from the array.
In [71]:
v2 = m2[:,2]
v2[2] = 999.
Out[71]:
In [72]:
m2
Out[72]:
In [73]:
v2
Out[73]:
The sub function produces an array view which accesses the actual elements in the original array.
In [74]:
v2 = sub(m2,:,2)
Out[74]:
In [75]:
v2[2] = 1.
Out[75]:
In [76]:
m2
Out[76]:
The sub function is preferred over explicit indexing if you need a subarray that you are not going to modify (sub does not create unnecessary copies) and for when you really do want to change the elements of a subarray.
Most linear algebra in Julia is performed on Matrix{Float64} objects using the LAPACK and BLAS software. By default Julia is compiled against OpenBLAS, an accelerated multi-threaded BLAS that also provides accelerated versions of some of the LAPACK functions. For example the LU decomposition and solutions of linear systems using LU are accelerated because they constitute the LAPACK benchmark. Julia can also be compiled with Intel's MKL (Math Kernel Library) BLAS, when they are available, to give accelerated performance. Understandably these are tuned for best performance on Intel processors.
Julia provides thin wrappers around many LAPACK and BLAS routines if you want maximal control. These mimic the BLAS and LAPACK routines without the annoyingly redundant arguments. Note that many internal linear algebra functions provide mutating versions.
In [77]:
whos(BLAS)
Operations like A'B and A\b call functions with names like Ac_mul_B and A_ldiv_B. There are mutating versions of these functions. The Ac denotes A' which, in general, is the conjugate transpose of A. For complex A there is a distinction between A', the conjugate transpose, and A.', the transpose. Usually you want A'. For real matrices A' and A.' are the same.
Notice that the type of cf33[:U] from the first Cholesky example is
In [78]:
typeof(cholfact(m33)[:U])
Out[78]:
Many operations can be steamlined if it is known that one of the operands is a triangular matrix or, even more specific, a unit triangular matrices. There are also special types for Symmetric and Diagonal matrices. Ongoing work will incorporate these matrix types more deeply into the Julia system.
When fitting a linear mixed model, the core computation can be expressed as solving a penalized least squares problem. In the MixedModels package there are several different PLSSolver types defined to facilitate this computation. When there is only a single random-effects term the system to be solved is block diagonal.
The PLSOne type is defined as
type PLSOne <: PLSSolver # Solver for models with a single random-effects term
Ad::Array{Float64,3} # diagonal blocks
Ab::Array{Float64,3} # base blocks
At::Symmetric{Float64} # lower right block
Ld::Array{Float64,3} # diagonal blocks
Lb::Array{Float64,3} # base blocks
Lt::Base.LinAlg.Cholesky{Float64} # lower right triangle
Zt::SparseMatrixCSC
end
and is updated from a triangular matrix λ as
## update!(s,lambda)->s : update Ld, Lb and Lt
function update!(s::PLSOne,λ::Triangular)
m,n,l = size(s)
n == size(λ,1) || error("Dimension mismatch")
Lt = copy!(s.Lt.UL,s.At.S)
Lb = copy!(s.Lb,s.Ab)
if n == 1 # shortcut for 1×1 λ
lam = λ[1,1]
Ld = map!(x -> sqrt(x*lam*lam + 1.), s.Ld, s.Ad)
Lb = scale!(reshape(Lb,(m,n*l)),lam ./ vec(Ld))
BLAS.syrk!('L','N',-1.0,Lb,1.0,Lt)
else
Ac_mul_B!(λ,reshape(copy!(s.Ld,s.Ad),(n,n*l)))
for k in 1:l
wL = A_mul_B!(sub(s.Ld,:,:,k),λ)
for j in 1:n # Inflate the diagonal
wL[j,j] += 1.
end
_, info = LAPACK.potrf!('L',wL) # i'th diagonal block of L
info == 0 || error("Cholesky failure at L diagonal block $k")
Base.LinAlg.A_rdiv_Bc!(A_mul_B!(sub(s.Lb,:,:,k),λ),Triangular(wL,:L,false))
end
BLAS.syrk!('L','N',-1.0,reshape(Lb,(m,n*l)),1.,Lt)
end
_, info = LAPACK.potrf!('L',Lt)
info == 0 || error("downdated X'X is not positive definite")
s
end
Sparse matrices and sparse vectors are available in the base Julia system but their description will need to wait for another time.
In [79]:
?sparse
Out[79]: