Merge pull request #3 from jiahao/master

merge back master

Author: dlfivefifty

.travis.yml

@@ -1,21 +1,13 @@
-language: cpp
-compiler:
-    - clang
+language: julia
+os:
+    - osx
+    - linux
+julia:
+    - release
+    - nightly
 notifications:
     email: false
-env:
-    - JULIAVERSION="juliareleases" 
-    - JULIAVERSION="julianightlies" 
-before_install:
-    - sudo add-apt-repository ppa:staticfloat/julia-deps -y
-    - sudo add-apt-repository ppa:staticfloat/${JULIAVERSION} -y
-    - sudo apt-get update -qq -y
-    - sudo apt-get install libpcre3-dev julia -y
-    - git config --global user.name "Travis User"
-    - git config --global user.email "travis@example.net"
-    - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
 script:
-    - julia -e 'Pkg.init(); Pkg.clone(pwd())'
-    - julia -e 'Pkg.test("RandomMatrices", coverage=true)'
-after_success:
-    - julia -e 'cd(Pkg.dir("RandomMatrices")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
+    - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
+    - julia -e 'Pkg.clone(pwd()); Pkg.build("RandomMatrices"); Pkg.test("RandomMatrices"; coverage=true)';
+    - julia -e 'cd(Pkg.dir("RandomMatrices")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())';

README.md

@@ -3,6 +3,11 @@ RandomMatrices.jl
 
 Random matrix package for [Julia](http://julialang.org).
 
+[![RandomMatrices on Julia release](http://pkg.julialang.org/badges/RandomMatrices_release.svg)](http://pkg.julialang.org/?pkg=RandomMatrices&ver=release)
+[![RandomMatrices on Julia nightly](http://pkg.julialang.org/badges/RandomMatrices_nightly.svg)](http://pkg.julialang.org/?pkg=RandomMatrices&ver=nightly)
+[![Build Status](https://travis-ci.org/jiahao/RandomMatrices.jl.png?branch=master)](https://travis-ci.org/jiahao/RandomMatrices.jl)
+[![Coverage Status](https://coveralls.io/repos/jiahao/RandomMatrices.jl/badge.svg?branch=master)](https://coveralls.io/r/jiahao/RandomMatrices.jl?branch=master)
+
 This extends the [Distributions](https://github.com/JuliaStats/Distributions.jl)
 package to provide methods for working with matrix-valued random variables,
 a.k.a. random matrices. State of the art methods for computing random matrix
@@ -11,11 +16,6 @@ samples and their associated distributions are provided.
 The names of the various ensembles can vary widely across disciplines. Where possible,
 synonyms are listed.
 
-[![Build Status](https://travis-ci.org/jiahao/RandomMatrices.jl.png?branch=master)](https://travis-ci.org/jiahao/RandomMatrices.jl)
-[![RandomMatrices](http://pkg.julialang.org/badges/RandomMatrices_release.svg)](http://pkg.julialang.org/?pkg=RandomMatrices&ver=release)
-[![RandomMatrices](http://pkg.julialang.org/badges/RandomMatrices_nightly.svg)](http://pkg.julialang.org/?pkg=RandomMatrices&ver=nightly)
-[![Coverage Status](https://img.shields.io/coveralls/jiahao/RandomMatrices.jl.svg)](https://img.shields.io/coveralls/jiahao/RandomMatrices.jl.svg)
-
 Additional functionality is provided when these optional packages are installed:
 - Symbolic manipulation of Haar matrices with [GSL.jl](https://github.com/jiahao/GSL.jl)
 - Invariant ensembles with [ApproxFun.jl](https://github.com/dlfivefifty/ApproxFun.jl)
@@ -71,7 +71,7 @@ Hermite, Laguerre(m) and Jacobi(m1, m2) ensembles.
   `GaussianLaguerreTridiagonalMatrix(n, m, beta)`,
   `GaussianJacobiSparseMatrix(n, m1, m2, beta)`
   each construct a sparse `n`x`n` matrix for the corresponding matrix ensemble
-  for arbitrary positive finite `beta`. 
+  for arbitrary positive finite `beta`.
   `GaussianHermiteTridiagonalMatrix(n, Inf)` is also allowed.
   These sampled matrices have the same eigenvalues as above but are much faster
   to diagonalize oweing to their sparsity. They also extend Dyson's threefold
@@ -87,7 +87,7 @@ Hermite, Laguerre(m) and Jacobi(m1, m2) ensembles.
    is applied to the raw QR decomposition. By default, `correction=1` (Edelman's correction) is
    used. Other valid values are `0` (no correction) and `2` (Mezzadri's correction).
  - `NeedsPiecewiseCorrection()` implements a simple test to see if a correction is necessary.
- 
+
 - `InvariantEnsemble(str,n)`
    Generates a unitary invariant ensemble, where str determines the
    potential of the ensemble, see below.
@@ -151,15 +151,15 @@ Provides finite-dimensional matrix representations of stochastic operators.
 In the following, `dt` is the time interval being discretized over and `t_end` is the final time.
 
 - `BrownianProcess(dt, t_end)` generates a vector corresponding to a Brownian random walk starting
-   from time `t=0` and position `x=0` 
-- `WhiteNoiceProcess(dt, t_end)` generates a vector corresponding to white noise.
+   from time `t=0` and position `x=0`
+- `WhiteNoiseProcess(dt, t_end)` generates a vector corresponding to white noise.
 - `StochasticAiryProcess(dt, t_end, beta)` generates the largest eigenvalue corresponding to the
    stochastic Airy process with real positive `beta`. This is known to be distributed in the `t_end -> Inf`
    limit to the `beta`-Tracy-Widom law.
-   
+
 # Invariant ensembles
 
-`InvariantEnsemble(str,n)` supports n x n unitary invariant ensemble 
+`InvariantEnsemble(str,n)` supports n x n unitary invariant ensemble
  with distribution
 
 `exp(- Tr Q(M)) dM`

REQUIRE

@@ -1,4 +1,5 @@
 julia 0.3
 Combinatorics
+Docile
 Distributions
 ODE

demos/book/8/randomgrowth.jl

@@ -7,15 +7,15 @@ function random_growth(M, N, q)
 
     G[1,1] = true
     # update the possible sets
-    imagesc((0,N), (M,0), G)
+    display(imagesc((0,N), (M,0), G))
     for t = 1:T
         sets = next_possible_squares(G)
         ## Highlights all the possible squares
         for i = 1:length(sets)
             idx = sets[i]::(Int,Int)
             G[idx[1], idx[2]] = 0.25
         end
-        imagesc((0,N), (M,0), G)
+        display(imagesc((0,N), (M,0), G))
         sleep(.01)
 
         ## Actual growth
@@ -26,12 +26,12 @@ function random_growth(M, N, q)
                 G[idx[1], idx[2]] = ison
             end
         end
-        imagesc((0,N), (M,0), G)
+        display(imagesc((0,N), (M,0), G))
         G[G .== 0.5] = 1
         G[G .== 0.25] = 0
         sleep(.01)
     end
-    G
+    return G
 end
 
 function next_possible_squares(G)

demos/book/8/randomgrowth2.jl

@@ -1,13 +1,14 @@
 #randomgrowth2.jl
-using Winston
+using PyPlot
 using ODE
-include("../1/gradient.jl")
 
 function randomgrowth2()
     num_trials = 2000
     g = 1
     q = 0.7
-    Ns = 80 # [50, 100, 200]
+    Ns = 80
+    # Ns = [50, 100, 200]
+    clf()
     # [-1.771086807411  08131947928329]
     # Gap = c d^N
     for jj = 1:length(Ns)
@@ -23,24 +24,25 @@ function randomgrowth2()
         end
         C = (B - N*om(g,q)) / (sigma_growth(g,q)*N^(1/3))
         d = 0.2
-        x, n = hist(C - exp(-N*0.0025 - 3), -6:d:4)
-        p = FramedPlot()
-        h = Histogram(n/(d*num_trials), step(x))
-        h.x0 = first(x)
-        add(p, h)
+        subplot(1,length(Ns),jj)
+        plt.hist(C - exp(-N*0.0025 - 3), normed = true)
+        xlim(-6, 4)
 
         ## Theory
         t0 = 4
         tn = -6
         dx = 0.005
         deq = (t, y) -> [y[2]; t*y[1]+2*y[1]^3; y[4]; y[1]^2]
         y0 = [airy(t0); airy(1,t0); 0; airy(t0)^2]      # boundary conditions
-        t, y = ode23(deq, t0:-dx:tn, y0)                # solve
-        F2 = exp(-y[:,3])                               # the distribution
+        t, y = ode23(deq, y0, t0:-dx:tn)                # solve
+        F2 = Float64[exp(-y[i][3]) for i = 1:length(y)]                               # the distribution
         f2 = gradient(F2, t)                            # the density
-       
-        add(p, Curve(t, f2, "color", "red", "linewidth", 3))
-        Winston.display(p)
+
+        # add(p, Curve(t, f2, "color", "red", "linewidth", 3))
+        # Winston.display(p)
+        subplot(1,length(Ns),jj)
+        plot(t, f2, "r", linewidth = 3)
+        ylim(0, 0.6)
         println(mean(C))
     end
 end
@@ -49,7 +51,7 @@ function G(N, M, q)
     # computes matrix G[N,M] for a given q
     GG = floor(log(rand(N,M))/log(q))
     #float(rand(N,M) < q)
-    
+
     # Compute the edges: the boundary
     for ii = 2:N
         GG[ii, 1] = GG[ii, 1] + GG[ii-1, 1]
@@ -64,7 +66,7 @@ function G(N, M, q)
     ## Plot
     # imagesc((0,N),(M,0),GG)
     # sleep(.01)
-    
+
     return  GG[N, M]
 end
 

src/HaarMeasure.jl

@@ -70,3 +70,44 @@ end
 #Zyczkowski and Kus, Random unitary matrices, J. Phys. A: Math. Gen. 27,
 #4235–4245 (1994).
 
+### Stewarts algorithm for n^2 orthogonal random matrix
+using Base.LinAlg: BlasInt
+for (s, elty) in (("dlarfg_", Float64),
+                  ("zlarfg_", Complex128))
+    @eval begin
+        function larfg!(n::Int, α::Ptr{$elty}, x::Ptr{$elty}, incx::Int, τ::Ptr{$elty})
+            ccall(($(Base.blasfunc(s)), Base.liblapack_name), Void,
+                (Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}),
+                &n, α, x, &incx, τ)
+        end
+    end
+end
+
+function Stewart(::Type{Float64}, n)
+    τ = Array(Float64, n)
+    H = randn(n, n)
+
+    pτ = pointer(τ)
+    pβ = pointer(H)
+    pH = pointer(H, 2)
+
+    for i = 0:n-2
+        larfg!(n - i, pβ + (n + 1)*8i, pH + (n + 1)*8i, 1, pτ + 8i)
+    end
+    LinAlg.QRPackedQ(H,τ)
+end
+function Stewart(::Type{Complex128}, n)
+    τ = Array(Complex128, n)
+    H = complex(randn(n, n), randn(n,n))
+
+    pτ = pointer(τ)
+    pβ = pointer(H)
+    pH = pointer(H, 2)
+
+    for i = 0:n-2
+        larfg!(n - i, pβ + (n + 1)*16i, pH + (n + 1)*16i, 1, pτ + 16i)
+    end
+    LinAlg.QRPackedQ(H,τ)
+end
+
+

src/RandomMatrices.jl

@@ -2,6 +2,10 @@ module RandomMatrices
 importall Distributions
 using Combinatorics
 
+if VERSION < v"0.4-"
+    using Docile
+end
+
 import Base.rand
 
 #If the GNU Scientific Library is present, turn on additional functionality.
@@ -32,6 +36,9 @@ include("densities/TracyWidom.jl")
 # Ginibre
 include("Ginibre.jl")
 
+# determinantal point processes
+include("dpp.jl")
+
 #Generating matrices of Haar measure
 include("Haar.jl")
 include("HaarMeasure.jl")

src/dpp.jl

@@ -0,0 +1,87 @@
+# Random samples from determinantal point processes
+
+@doc doc"""
+Computes a random sample from the determinantal point process defined by the
+spectral factorization object `L`.
+
+Inputs:
+
+    `L`: `Eigen` factorization object of an N x N matrix
+
+Output:
+
+    `Y`: A `Vector{Int}` with entries in [1:N].
+
+References:
+
+    Algorithm 18 of \cite{HKPV05}, as described in Algorithm 1 of \cite{KT12}.
+
+    @article{HKPV05,
+        author = {Hough, J Ben and Krishnapur, Manjunath and Peres, Yuval and Vir\'{a}g, B\'{a}lint},
+        doi = {10.1214/154957806000000078},
+        journal = {Probability Surveys},
+        pages = {206--229},
+        title = {Determinantal Processes and Independence},
+        volume = {3},
+        year = {2005}
+        archivePrefix = {arXiv},
+        eprint = {0503110},
+    }
+
+    @article{KT12,
+        author = {Kulesza, Alex and Taskar, Ben},
+        doi = {10.1561/2200000044},
+        journal = {Foundations and Trends in Machine Learning},
+        number = {2-3},
+        pages = {123--286},
+        title = {Determinantal Point Processes for Machine Learning},
+        volume = {5},
+        year = {2012},
+        archivePrefix = {arXiv},
+        eprint = {1207.6083},
+    }
+""" ->
+function rand{S<:Real,T}(L::Base.LinAlg.Eigen{S,T})
+    N = length(L.values)
+    J = Int[]
+    for n=1:N
+        λ = L.values[n]
+        rand() < λ/(λ+1) && push!(J, n)
+    end
+
+    V = L.vectors[:, J]
+    Y = Int[]
+    nV = size(V, 2)
+    while true
+        # Select i from 𝒴=[1:N] (ground set) with probabilities
+        # Pr(i) = 1/|V| Σ_{v∈V} (v⋅eᵢ)²
+
+        #Compute selection probabilities
+        Pr = zeros(N)
+        for i=1:N
+            for j=1:nV #TODO this loop is a bottleneck - why?
+                Pr[i] += (V[i,j])^2 #ith entry of jth eigenvector
+            end
+            Pr[i] /= nV
+        end
+        @assert abs(1-sum(Pr)) < N*eps() #Check normalization
+
+        #Simple discrete sampler
+        i, ρ = N, rand()
+        for j=1:N
+            if ρ < Pr[j]
+                i = j
+                break
+            else
+                ρ -= Pr[j]
+            end
+        end
+        push!(Y, i)
+        nV == 1 && break #Done
+        #V = V⊥ #an orthonormal basis for the subspace of V ⊥ eᵢ
+        V[i, :] = 0 #Project out eᵢ
+        V = full(qrfact!(V)[:Q])[:, 1:nV-1]
+        nV = size(V, 2)
+    end
+    Y
+end

test/Haar.jl

@@ -12,4 +12,9 @@ Q=UniformHaar(2, N)
 println("Case 3")
 println("E(A*Q*B*Q'*A*Q*B*Q') = ", eval(expectation(N, :Q, :(A*Q*B*Q'*A*Q*B*Q'))))
 
+for elty in (Float64, Complex128)
+	A = Stewart(elty, N)
+	@test_approx_eq A'A eye(N)
+end
+
 end #_HAVE_GSL