{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MatrixPerspective.jl" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Julia Version 1.5.1\n", "Commit 697e782ab8 (2020-08-25 20:08 UTC)\n", "Platform Info:\n", " OS: macOS (x86_64-apple-darwin19.5.0)\n", " CPU: Intel(R) Core(TM) i5-8279U CPU @ 2.40GHz\n", " WORD_SIZE: 64\n", " LIBM: libopenlibm\n", " LLVM: libLLVM-9.0.1 (ORCJIT, skylake)\n" ] } ], "source": [ "versioninfo()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Timing\n", "\n", "The following code compares the three methods, interior point method (by using MOSEK, for $p < 100$), bisection, and semismooth Newton, for finding the unique root of the semismooth function\n", "\n", "$f(\\mu) = 1 - \\mathbf{e}^T(\\bar{\\mathbf{X}} - \\mu \\mathbf{e}\\mathbf{e}^T)_{+}\\mathbf{e}$\n", "\n", "where $\\mathbf{e}=(0, \\dotsc, 0, 1)^T$ and\n", "\n", "$\\bar{\\mathbf{X}} = \\begin{bmatrix} -\\mathbf{X} & \\frac{1}{\\sqrt{2}}\\mathbf{y} \\\\\n", " \\frac{1}{\\sqrt{2}}\\mathbf{y}^T & 1 \\end{bmatrix}$\n", " \n", "for given input $(\\mathbf{X}, \\mathbf{y})$. From this root the prox operator\n", "\n", "$\\mathrm{prox}_{\\phi}(\\mathbf{X}, \\mathbf{y}),\n", " \\quad\n", " \\phi(\\boldsymbol{\\Omega}, \\boldsymbol{\\eta}) = \\begin{cases}\n", " \\frac{1}{2}\\boldsymbol{\\eta}^T\\boldsymbol{\\Omega}^{\\dagger}\\boldsymbol{\\eta}, &\n", " \\boldsymbol{\\Omega} \\succeq \\mathbf{0},~\\boldsymbol{\\eta} \\in \\mathcal{R}(\\boldsymbol{\\Omega}) \\\\\n", " \\infty, & \\text{otherwise}\n", " \\end{cases}$\n", "\n", "can be computed in a closed form." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first table reports the mean of the performance measures, and the second table contains the standard deviation. The $p=5$ case (especially for MOSEK) should be ignored since there is an overhead of JIT compilation of the code." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "using MatrixPerspective\n", "\n", "using LinearAlgebra\n", "import LinearAlgebra.BLAS.BlasInt\n", "\n", "using Random\n", "using DataFrames, Statistics\n", "\n", "Random.seed!(1234);\n", "\n", "reps = 10 #reps = 100\n", "tol = 1e-8 # default\n", "\n", "# KKT measures |g'(mu)|\n", "df = DataFrame(p=Int[], Method=String[], Iters=Float64[], Secs=Float64[], KKT=Float64[], Obj=Float64[])\n", "# MOSEK stalls if p > 30\n", "dims = [5, 10, 30, 50, 100, 500, 1000, 2000] \n", "#dims = [10, 30] #dims = [100] #dims = [1000]\n", "nummethods = 2\n", "Means = zeros(nummethods, length(dims), size(df)[2]);\n", "mosek_maxdim = 100\n", "for i = 1:length(dims)\n", "\tp = dims[i]\n", "\tn = p + 2\n", "\n", "\te = zeros(p + 1) # last elementary unit vector\n", "\te[end] = 1\n", "\n", "\t# workspaces\n", "\tQ = Matrix{Float64}(undef, n, n)\n", "\tevec = Vector{Float64}(undef, n)\n", "\n", "\td = Vector{Float64}(undef, n)\n", "\n", "\tindxq = Vector{BlasInt}(undef, n)\n", "\n", "\tz = Vector{Float64}(undef, n)\n", "\tdlambda = Vector{Float64}(undef, n)\n", "\tw = Vector{Float64}(undef, n)\n", "\tQ2 = Vector{Float64}(undef, n * n)\n", "\n", "\tindx = Vector{BlasInt}(undef, n + 1)\n", "\tindxc = Vector{BlasInt}(undef, n)\n", "\tindxp = Vector{BlasInt}(undef, n)\n", "\tcoltyp = Vector{BlasInt}(undef, n)\n", "\n", "\tS = Vector{Float64}(undef, n * n )\n", "\n", "\tM = zeros(n - 1, n - 1) # for second derivative\n", "\n", "\talph = Vector{Bool}(undef, n - 1)\n", "\tgam = similar(alph)\n", "\tbet = similar(alph)\n", "\n", "println(\"p = \", p)\n", "\tfor rep=1:reps\n", "\t\tX = Matrix(Symmetric(randn(p, p)))\n", "\t\ty = randn(p)\n", "\t\twhile !(minimum(eigvals(X + 0.5 * y * y')) < 0)\n", "\t\t\tX = Diagonal(randn(p))\n", "\t\t\ty = randn(p)\n", "\t\tend\n", "\t\tG = [X -y/sqrt(2); -y'/sqrt(2) 1]\n", "\n", "\t\tif p <= mosek_maxdim\n", "\t\t\tsecs = @elapsed optval, primal, dual = mosek_ls_cov_adj(Matrix(G))\n", "\t\t\tmu = dual[2].dual\n", "\n", "\t\t\tlam2, P2 = eigen(G + mu * e * e')\n", "\t\t\tLam2 = P2 * Diagonal(max.(lam2, 0)) * P2'\n", "\t\t\tmosek_deriv1 = 1 - Lam2[end, end]\n", "\t\telse\n", "\t\t\t(optval, mu) = (NaN, NaN) # if MOSEK stalls\n", "\t\t\tmosek_deriv1 = NaN\n", "\t\t\tsecs = NaN\n", "\t\tend\n", "\t\tpush!(df, [p, \"MOSEK\", NaN, secs, abs(mosek_deriv1), optval])\n", "\n", "\t\tsecs = @elapsed nu, C, convhist = bisect!(Matrix(G), Q, evec, d,\n", "\t\t\t\t\t\t\t\t\t\tindxq, z, dlambda, w, Q2,\n", "\t\t\t\t\t\t\t\t\t\tindx, indxc, indxp, coltyp, S,\n", "\t\t\t\t\t\t\t\t\t\tM, alph, bet, gam; maxiter=50)\n", "\n", "\t\t# reconstruct the dual optimal value\n", "\t\tdualval = -0.5 * norm(C, 2)^2 + nu + 0.5 * norm(G, 2)^2\n", "\n", "\t\tderiv1 = 1 - C[end, end]\n", "\n", "\t\tpush!(df, [p, \"Bisection\", convhist.iters, secs, abs(deriv1), dualval])\n", "\n", "\t\tsecs = @elapsed nu, C, convhist = dual_ls_cov_adj!(Matrix(G), Q, evec, d,\n", "\t\t\t\t\t\t\t\t\t\tindxq, z, dlambda, w, Q2,\n", "\t\t\t\t\t\t\t\t\t\tindx, indxc, indxp, coltyp, S,\n", "\t\t\t\t\t\t\t\t\t\tM, alph, bet, gam; maxiter=50)\n", "\n", "\t\t# reconstruct the dual optimal value\n", "\t\tdualval = -0.5 * norm(C, 2)^2 + nu + 0.5 * norm(G, 2)^2\n", "\n", "\t\tderiv1 = 1 - C[end, end]\n", "\n", "\t\tpush!(df, [p, \"Newton\", convhist.iters, secs, abs(deriv1), dualval])\n", "\tend\n", "end\n", "gdf = groupby(df, [:p, :Method])\n", "mdf = combine(gdf, names(gdf)[3:end] .=> mean)\n", "println(mdf)\n", "sdf = combine(gdf, names(gdf)[3:end] .=> std)\n", "println(sdf)\n", "\n" ] } ], "source": [ ";cat timing.jl" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p = 5\n", "p = 10\n", "p = 30\n", "p = 50\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "┌ Warning: Problem status SLOW_PROGRESS; solution may be inaccurate.\n", "└ @ Convex /Users/jhwon/.julia/packages/Convex/aYxJA/src/solution.jl:252\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "p = 100\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "┌ Warning: Problem status SLOW_PROGRESS; solution may be inaccurate.\n", "└ @ Convex /Users/jhwon/.julia/packages/Convex/aYxJA/src/solution.jl:252\n", "┌ Warning: Problem status SLOW_PROGRESS; solution may be inaccurate.\n", "└ @ Convex /Users/jhwon/.julia/packages/Convex/aYxJA/src/solution.jl:252\n", "┌ Warning: Problem status SLOW_PROGRESS; solution may be inaccurate.\n", "└ @ Convex /Users/jhwon/.julia/packages/Convex/aYxJA/src/solution.jl:252\n", "┌ Warning: Problem status SLOW_PROGRESS; solution may be inaccurate.\n", "└ @ Convex /Users/jhwon/.julia/packages/Convex/aYxJA/src/solution.jl:252\n", "┌ Warning: Problem status SLOW_PROGRESS; solution may be inaccurate.\n", "└ @ Convex /Users/jhwon/.julia/packages/Convex/aYxJA/src/solution.jl:252\n", "┌ Warning: Problem status SLOW_PROGRESS; solution may be inaccurate.\n", "└ @ Convex /Users/jhwon/.julia/packages/Convex/aYxJA/src/solution.jl:252\n", "┌ Warning: Problem status SLOW_PROGRESS; solution may be inaccurate.\n", "└ @ Convex /Users/jhwon/.julia/packages/Convex/aYxJA/src/solution.jl:252\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "p = 500\n", "p = 1000\n", "p = 2000\n", "24×6 DataFrame\n", "│ Row │ p │ Method │ Iters_mean │ Secs_mean │ KKT_mean │ Obj_mean │\n", "│ │ \u001b[90mInt64\u001b[39m │ \u001b[90mString\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │\n", "├─────┼───────┼───────────┼────────────┼─────────────┼─────────────┼───────────┤\n", "│ 1 │ 5 │ MOSEK │ NaN │ 2.49323 │ 1.21816e-5 │ 6.53939 │\n", "│ 2 │ 5 │ Bisection │ 28.4 │ 0.157635 │ 5.97528e-9 │ 6.53939 │\n", "│ 3 │ 5 │ Newton │ 4.2 │ 0.0943318 │ 5.58998e-12 │ 6.53939 │\n", "│ 4 │ 10 │ MOSEK │ NaN │ 0.00983312 │ 1.828e-5 │ 27.5118 │\n", "│ 5 │ 10 │ Bisection │ 29.4 │ 0.000295918 │ 4.56144e-9 │ 27.5118 │\n", "│ 6 │ 10 │ Newton │ 4.9 │ 0.000165922 │ 2.34184e-10 │ 27.5118 │\n", "│ 7 │ 30 │ MOSEK │ NaN │ 0.102232 │ 4.92603e-6 │ 225.301 │\n", "│ 8 │ 30 │ Bisection │ 31.8 │ 0.001258 │ 5.44707e-9 │ 225.301 │\n", "│ 9 │ 30 │ Newton │ 5.7 │ 0.000554437 │ 1.18062e-9 │ 225.301 │\n", "│ 10 │ 50 │ MOSEK │ NaN │ 0.657758 │ 7.47059e-6 │ 654.111 │\n", "│ 11 │ 50 │ Bisection │ 33.2 │ 0.00311099 │ 4.26987e-9 │ 654.111 │\n", "│ 12 │ 50 │ Newton │ 5.8 │ 0.00112343 │ 5.79867e-10 │ 654.111 │\n", "│ 13 │ 100 │ MOSEK │ NaN │ 21.2737 │ 6.95149e-6 │ 2496.79 │\n", "│ 14 │ 100 │ Bisection │ 34.5 │ 0.0123491 │ 5.21449e-9 │ 2496.79 │\n", "│ 15 │ 100 │ Newton │ 6.5 │ 0.00433008 │ 2.44491e-9 │ 2496.79 │\n", "│ 16 │ 500 │ MOSEK │ NaN │ NaN │ NaN │ NaN │\n", "│ 17 │ 500 │ Bisection │ 36.9 │ 0.331822 │ 4.0035e-9 │ 62564.1 │\n", "│ 18 │ 500 │ Newton │ 7.9 │ 0.109068 │ 9.81724e-10 │ 62564.1 │\n", "│ 19 │ 1000 │ MOSEK │ NaN │ NaN │ NaN │ NaN │\n", "│ 20 │ 1000 │ Bisection │ 37.8 │ 2.19949 │ 4.43912e-9 │ 2.50248e5 │\n", "│ 21 │ 1000 │ Newton │ 8.0 │ 0.740104 │ 3.05567e-12 │ 2.50248e5 │\n", "│ 22 │ 2000 │ MOSEK │ NaN │ NaN │ NaN │ NaN │\n", "│ 23 │ 2000 │ Bisection │ 39.6 │ 14.6615 │ 5.84813e-9 │ 1.00103e6 │\n", "│ 24 │ 2000 │ Newton │ 8.0 │ 4.32247 │ 2.88162e-9 │ 1.00103e6 │\n", "24×6 DataFrame\n", "│ Row │ p │ Method │ Iters_std │ Secs_std │ KKT_std │ Obj_std │\n", "│ │ \u001b[90mInt64\u001b[39m │ \u001b[90mString\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │\n", "├─────┼───────┼───────────┼───────────┼─────────────┼─────────────┼─────────┤\n", "│ 1 │ 5 │ MOSEK │ NaN │ 7.8583 │ 1.28047e-5 │ 3.22705 │\n", "│ 2 │ 5 │ Bisection │ 0.966092 │ 0.498075 │ 3.69496e-9 │ 3.22705 │\n", "│ 3 │ 5 │ Newton │ 0.632456 │ 0.298139 │ 9.18807e-12 │ 3.22705 │\n", "│ 4 │ 10 │ MOSEK │ NaN │ 0.000569673 │ 1.72334e-5 │ 5.87331 │\n", "│ 5 │ 10 │ Bisection │ 1.95505 │ 3.38437e-5 │ 3.18761e-9 │ 5.87331 │\n", "│ 6 │ 10 │ Newton │ 0.316228 │ 4.13268e-5 │ 4.37243e-10 │ 5.87331 │\n", "│ 7 │ 30 │ MOSEK │ NaN │ 0.012398 │ 7.51272e-6 │ 25.0831 │\n", "│ 8 │ 30 │ Bisection │ 1.0328 │ 4.06005e-5 │ 3.51572e-9 │ 25.0831 │\n", "│ 9 │ 30 │ Newton │ 0.483046 │ 2.50481e-5 │ 3.14195e-9 │ 25.0831 │\n", "│ 10 │ 50 │ MOSEK │ NaN │ 0.110312 │ 9.48564e-6 │ 28.42 │\n", "│ 11 │ 50 │ Bisection │ 0.918937 │ 8.7696e-5 │ 2.88447e-9 │ 28.42 │\n", "│ 12 │ 50 │ Newton │ 0.421637 │ 6.49869e-5 │ 1.49951e-9 │ 28.42 │\n", "│ 13 │ 100 │ MOSEK │ NaN │ 4.43358 │ 6.50369e-6 │ 60.6632 │\n", "│ 14 │ 100 │ Bisection │ 0.707107 │ 0.000243857 │ 3.36907e-9 │ 60.6632 │\n", "│ 15 │ 100 │ Newton │ 0.527046 │ 0.000303794 │ 3.5232e-9 │ 60.6632 │\n", "│ 16 │ 500 │ MOSEK │ NaN │ NaN │ NaN │ NaN │\n", "│ 17 │ 500 │ Bisection │ 0.875595 │ 0.00846495 │ 2.54059e-9 │ 436.492 │\n", "│ 18 │ 500 │ Newton │ 0.316228 │ 0.00451969 │ 3.10448e-9 │ 436.492 │\n", "│ 19 │ 1000 │ MOSEK │ NaN │ NaN │ NaN │ NaN │\n", "│ 20 │ 1000 │ Bisection │ 1.61933 │ 0.0865994 │ 3.93659e-9 │ 618.872 │\n", "│ 21 │ 1000 │ Newton │ 0.0 │ 0.0137851 │ 2.3702e-12 │ 618.872 │\n", "│ 22 │ 2000 │ MOSEK │ NaN │ NaN │ NaN │ NaN │\n", "│ 23 │ 2000 │ Bisection │ 0.966092 │ 0.378565 │ 3.164e-9 │ 1197.17 │\n", "│ 24 │ 2000 │ Newton │ 0.0 │ 0.152881 │ 5.28822e-10 │ 1197.17 │\n" ] } ], "source": [ "include(\"timing.jl\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## PDHG\n", "The following code illustrates how to use `prox_matrixperspective!()` for the PDHG algorithm for Gaussian joint likelihood estimation." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "using MatrixPerspective\n", "\n", "using LinearAlgebra\n", "import LinearAlgebra.BLAS.BlasInt\n", "\n", "using IterativeSolvers\n", "\n", "# Joint MLE of multivariate Gaussian natural parameters\n", "# PDHG code based on http://proximity-operator.net/tutorial.html\n", "\n", "function gaussianmle(X::Matrix{T}, cvec::Array{Vector{T},1}; \n", "\t\t\t\t\t ep::T = 10 / size(X)[2]^2, \n", "\t\t\t\t\t init::Matrix{T} = Array{Float64}(undef, 0, 0), \n", "\t\t\t\t\t maxiter::Integer = 1000, \n", "\t\t\t\t\t opttol::T = 1e-4, \n", "\t\t\t\t\t log::Bool = false\n", "\t\t\t\t\t) where T <: Float64\n", "\tn, p = size(X) # sample size, dimension\n", "\tm = length(cvec)\n", "\n", "\t# workspaces\n", "\t_n = p + 2\n", "\t_Q = Matrix{Float64}(undef, _n, _n)\n", "\t_evec = Vector{Float64}(undef, _n)\n", "\n", "\t_d = Vector{Float64}(undef, _n)\n", "\n", "\t_indxq = Vector{BlasInt}(undef, _n)\n", "\n", "\t_z = Vector{Float64}(undef, _n)\n", "\t_dlambda = Vector{Float64}(undef, _n)\n", "\t_w = Vector{Float64}(undef, _n)\n", "\t_Q2 = Vector{Float64}(undef, _n * _n)\n", "\n", "\t_indx = Vector{BlasInt}(undef, _n + 1)\n", "\t_indxc = Vector{BlasInt}(undef, _n)\n", "\t_indxp = Vector{BlasInt}(undef, _n)\n", "\t_coltyp = Vector{BlasInt}(undef, _n)\n", "\n", "\t_S = Vector{Float64}(undef, _n * _n )\n", "\n", "\t_M = zeros(_n - 1, _n - 1) # for second derivative\n", "\n", "\t_alph = Vector{Bool}(undef, _n - 1)\n", "\t_gam = similar(_alph)\n", "\t_bet = similar(_alph)\n", "\n", "\tsymmetrize!(A) = begin\n", " \tfor c in 1:size(A, 2)\n", " \t\tfor r in c:size(A, 1)\n", " \t\t\tA[r, c] += A[c, r]\n", "\t\t\t\tA[r, c] *= 0.5\n", " \t\t\tA[c, r] = A[r, c]\n", " \t\tend\n", " \tend\n", " \tA\n", "\tend\n", "\n", "\t# select the step-sizes\n", "\tL_f = 0.0 \t\t\t # Lipschitz modulus of $f$\n", "\tL_K = m \t\t\t\t# |K'K|_2\n", "\ttau = 2 / (L_f + 2) \n", "\tsigma = (1/tau - L_f/2) / L_K \n", "\n", "\tS = X' * X / n # second moment\n", "\tmu = reshape(sum(X, dims=1) / n, p) # sample mean\n", "\n", "\t# initialize the primal solution\n", "\tif isempty(init) \n", "\t\t#Omega = 1.0 * Matrix(I, p, p) # identity matrix\n", "\t\tOmega = Matrix(Symmetric(inv(S - mu * mu' + 1e-2 * I)))\n", "\telse \n", "\t\tOmega = init\n", "\tend\n", "\tOmega_old = similar(Omega)\n", "\t#Omega_old = Omega\n", "\t#eta = zeros(p)\n", "\teta = Omega * mu\n", "\t#eta_old = similar(eta)\n", "\teta_old = similar(eta)\n", "\t# initialize the dual solution\n", "\tY = [copy(Omega) for i=1:m+1]\n", "\tzeta = copy(eta)\n", "\n", "\t# execute the algorithm\n", "\thistory = ConvergenceHistory(partial = !log)\n", "\thistory[:tol] = opttol\n", "\tIterativeSolvers.reserve!(T, history, :objval, maxiter)\n", "\tIterativeSolvers.reserve!(Float64, history, :itertime, maxiter)\n", "\tIterativeSolvers.nextiter!(history)\n", "\tfor it = 1:maxiter\n", "\t\ttic = time()\n", "if it % 100 == 0 \n", "\tprintln(\"it = \", it)\n", "end\n", " \t# primal forward-backward step\n", "\t\tcopyto!(Omega_old, Omega)\n", "\t\tcopyto!(eta_old, eta)\n", "\t\t## prox opertor of -logdet\n", "\t\tM = (Omega - tau * (sum(Y) + S)) ./ ( 1 + ep * tau)\n", "\t\tOmD, OmP = eigen!(Symmetric(M))\n", "\t\tevals = 0.5 * (OmD + sqrt.(OmD.^2 .+ 4 * tau / (1 + ep * tau)))\n", "\t\tfill!(Omega, zero(T))\n", "\t\t@inbounds for k=1:length(evals)\n", " \t\t@views pvec = OmP[:, k]\n", " \t\tBLAS.gemm!('N', 'T', evals[k], pvec, pvec, 1.0, Omega)\n", " \tend\n", "\t\t## eta update\n", "\t\teta -= tau * (zeta - 2 * mu)\n", "\n", "\t\t# primal adjustment\n", "\t\tOmega_tilde = 2 * Omega - Omega_old\n", "\t\teta_tilde = 2 * eta - eta_old\n", " \n", " \t# dual forward-backward step\n", "\t\tfor i=1:m\n", "\t\t\t#Y_tilde = cvec[i] * cvec[i]' - Y[i] ./ sigma - Omega_tilde\n", "\t\t\t#getPSDpart!(Y_tilde)\n", "\t\t\t#Y[i] = -sigma * Y_tilde\n", "\t\t\t## less allocation using low-level matrix functions\n", "\t\t\t#Y[i] .+= sigma * Omega_tilde\n", "\t\t\tBLAS.axpy!(sigma, Omega_tilde, Y[i]) \n", "\t\t\t# Y[i] .= sigma * cvec[i] * cvex[i]' - Y[i]\n", "\t\t\tBLAS.gemm!('N', 'T', sigma, cvec[i], cvec[i], -1.0, Y[i])\n", "\t\t\t# ([Y[i])_+\n", "\t\t\tgetPSDpart!(Y[i])\n", "\t\t\t# Y[i] .*= -1.0\n", "\t\t\trmul!(Y[i], -1.0)\n", "\t\tend\n", "\t\t# prox operator of $\\phi^*$\n", "\t\tM = Y[m + 1] + sigma * Omega_tilde\n", "\t\tsymmetrize!(M)\n", "\t\tz = zeta + sigma * eta_tilde\n", "\t\t# the following puts *primal* variables to Y[m + 1] and zeta\n", "\t\tprox_matrixperspective!( Y[m + 1], zeta,\n", "\t\t\t\t\t\t\t\t\tM, z, 1.0,\n", "\t\t\t\t\t\t \t\t\t_Q, _evec, _d, \n", "\t\t\t\t\t\t\t\t\t_indxq, _z, _dlambda, _w, _Q2,\n", "\t\t\t\t\t\t\t\t\t_indx, _indxc, _indxp, _coltyp, _S,\n", "\t\t\t\t\t\t\t\t\t_M, _alph, _bet, _gam\n", "\t\t\t\t\t\t\t\t )\n", "\t\t## restore dual variables\n", "\t\t# Y[m + 1] .= M - Y[m + 1]\n", "\t\tBLAS.axpy!(-1.0, M, Y[m + 1]) # Y[m + 1] .= -M + Y[m + 1]\n", "\t\trmul!(Y[m + 1], -1.0) # Y[m + 1] .*= -1.0\n", "\t\t# zeta .= z - zeta\n", "\t\tBLAS.axpy!(-1.0, z, zeta) # zeta .= -z + zeta\n", "\t\trmul!(zeta, -1.0) # zeta .*= -1.0\n", "\n", " \t# objective value\n", "\t\tD, P = eigen(Symmetric(Omega))\n", "\t\tobjval = -sum(Base.log.(max.(D, 1e-15))) + tr(Omega * S) - 2 * dot(mu, eta) + (eta' * P) * (pinv(Diagonal(D)) * P' * eta) * 0.5 + 0.5 * ep * sum(abs2.(D))\n", "if it % 100 == 1\n", "\tprintln(\"objval = \", objval)\n", "end\n", "\t\ttoc = time()\n", "\t\tpush!(history, :objval, objval)\n", "\t\tpush!(history, :itertime, toc - tic)\n", " \t# stopping rule\n", " \tif norm(vec(Omega - Omega_old)) < opttol * norm(vec(Omega_old)) && \n", "\t\t norm(eta - eta_old) < opttol * norm(eta_old) && it > 10\n", "\t\t\tIterativeSolvers.setconv(history, true)\n", " \tbreak\n", " \tend\n", "\t\tIterativeSolvers.nextiter!(history)\n", "\tend # end for\n", "\tlog && IterativeSolvers.shrink!(history)\n", "\teta, Omega, Y, history\n", "end # end of function\n", "\n" ] } ], "source": [ ";cat gaussianmle.jl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Function `gaussianmle()` is called as follows." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "include(\"gaussianmle.jl\")\n", "using Random, LinearAlgebra, SparseArrays\n", "\n", "Random.seed!(123)\n", "n, p = 60, 100\n", "# data matrix\n", "pcov = 0.3 * ones(p, p) + 0.7 * I # covariance matrix\n", "pchol = cholesky(pcov)\n", "X = randn(n, p) * pchol.U # correlated predictors\n", "S = X' * X / n # second monent\n", "mu = reshape(sum(X, dims=1) / n , p ) # sample mean\n", "\n", "# Variance constraints\n", "# cvec should be scaled so that\n", "# cvec[i]' * (Omega \\ cvec[i]) <= 1\n", "cvec = Array{Vector{Float64}, 1}()\n", "m = 5\n", "for i=1:m\n", "\te = zeros(p)\n", "\te[i] = 1.0\n", "\tpush!(cvec, e)\n", "end\n", "\n", "maxit = 5000\n", "@time eta, Omega, Y, history = gaussianmle(X, cvec, maxiter=maxit, opttol=1e-5)\n", "\n" ] } ], "source": [ ";cat testgaussianmle.jl" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "objval = 824.1079588444118\n", "it = 100\n", "objval = -68.55345451618902\n", "it = 200\n", "objval = -90.7432608163429\n", "it = 300\n", "objval = -106.47412453177637\n", "it = 400\n", "objval = -119.48174383805825\n", "it = 500\n", "objval = -130.905231884441\n", "it = 600\n", "objval = -141.24124058414466\n", "it = 700\n", "objval = -150.75060287270733\n", "it = 800\n", "objval = -159.58819651238332\n", "it = 900\n", "objval = -167.8548291740341\n", "it = 1000\n", "objval = -175.62112227294827\n", "it = 1100\n", "objval = -182.93963829174476\n", "it = 1200\n", "objval = -189.85149972890522\n", "it = 1300\n", "objval = -196.39020970969847\n", "it = 1400\n", "objval = -202.58396017158287\n", "it = 1500\n", "objval = -208.45708168793593\n", "it = 1600\n", "objval = -214.03098642732363\n", "it = 1700\n", "objval = -219.32480171804784\n", "it = 1800\n", "objval = -224.3558092685453\n", "it = 1900\n", "objval = -229.13975914037536\n", "it = 2000\n", "objval = -233.6911010552277\n", "it = 2100\n", "objval = -238.02315987101295\n", "it = 2200\n", "objval = -242.1482724796969\n", "it = 2300\n", "objval = -246.07789742449864\n", "it = 2400\n", "objval = -249.82270476506264\n", "it = 2500\n", "objval = -253.39265129471005\n", "it = 2600\n", "objval = -256.7970446309472\n", "it = 2700\n", "objval = -260.0445986525942\n", "it = 2800\n", "objval = -263.14348205382316\n", "it = 2900\n", "objval = -266.1013613073703\n", "it = 3000\n", "objval = -268.92543899954114\n", "it = 3100\n", "objval = -271.62248826913253\n", "it = 3200\n", "objval = -274.1988839186122\n", "it = 3300\n", "objval = -276.660630647678\n", "it = 3400\n", "objval = -279.0133887725162\n", "it = 3500\n", "objval = -281.2624977290658\n", "it = 3600\n", "objval = -283.41299760919617\n", "it = 3700\n", "objval = -285.46964894025103\n", "it = 3800\n", "objval = -287.4369508881293\n", "it = 3900\n", "objval = -289.3191580396211\n", "it = 4000\n", "objval = -291.12029589988504\n", "it = 4100\n", "objval = -292.8441752243772\n", "it = 4200\n", "objval = -294.4944052907238\n", "it = 4300\n", "objval = -296.0744062042605\n", "it = 4400\n", "objval = -297.58742032087287\n", "it = 4500\n", "objval = -299.0365228620533\n", "it = 4600\n", "objval = -300.42463178951584\n", "it = 4700\n", "objval = -301.75451700000485\n", "it = 4800\n", "objval = -303.0288088951596\n", "it = 4900\n", "objval = -304.25000637601295\n", "it = 5000\n", "118.756086 seconds (13.64 M allocations: 17.783 GiB, 2.01% gc time)\n" ] }, { "data": { "text/plain": [ "([-116.11434098621383, 48.95213796995712, 90.65112808590104, 74.88200654160379, -109.65229977959252, -98.84490753993417, -30.09665095952194, 136.9517017393632, -104.99144451079381, -18.399580485564563 … -49.91108946571656, 117.59080136676128, -68.00851618779222, -94.5335199879541, 74.52083217154109, 90.99154106488278, 108.42836157637981, 133.7807753285948, -48.671229175687735, -102.01022091106664], [21.33766714582228 -2.6238047353834855 … 0.6781649941396397 4.81120439765941; -2.6238047353834855 12.248940051450315 … -2.527585193176585 -3.467548444667807; … ; 0.6781649941396397 -2.527585193176585 … 17.906439248985176 3.2836072748514855; 4.81120439765941 -3.467548444667807 … 3.2836072748514855 21.367264541224216], [[-0.14575876261675827 -0.03563277907757659 … -0.03404245270484068 -0.06075941502196815; -0.03563277907757659 -0.008710933888275197 … -0.008322156244423087 -0.014853493357741087; … ; -0.03404245270484068 -0.008322156244423087 … -0.007950730133517705 -0.014190567175007886; -0.06075941502196815 -0.014853493357741087 … -0.014190567175007886 -0.025327509972888063], [-0.008689324794789004 -0.03553955835073861 … -0.006459451495274717 -0.01291964301956979; -0.03553955835073861 -0.14535769321489891 … -0.026419320114234028 -0.05284166696589258; … ; -0.006459451495274717 -0.026419320114234028 … -0.00480181309885309 -0.009604176318880721; -0.01291964301956979 -0.05284166696589258 … -0.009604176318880721 -0.01920945294313533], [-0.0007944588520018045 -0.0017130282012015343 … -0.0008562326215157404 -0.0015634040994131277; -0.0017130282012015343 -0.003693665959813736 … -0.001846226049026201 -0.00337104345356667; … ; -0.0008562326215157404 -0.001846226049026201 … -0.0009228096588016267 -0.0016849678081576932; -0.0015634040994131277 -0.00337104345356667 … -0.0016849678081576932 -0.0030766003448800656], [-0.0 -0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0; … ; -0.0 -0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0], [-0.03234378615574314 -0.020965129193548145 … -0.011242839572479554 -0.026884814875212406; -0.020965129193548145 -0.013589523501846373 … -0.007287569334164546 -0.017426643080354694; … ; -0.011242839572479554 -0.007287569334164546 … -0.0039080595278443535 -0.00934527760981235; -0.026884814875212406 -0.017426643080354694 … -0.00934527760981235 -0.02234720658224793], [-0.010322316049599323 0.01708816716554684 … 0.00295235387583076 -0.008294524720988461; 0.01708816716554684 -0.028288753771398945 … -0.004887499696737629 0.013731242514756659; … ; 0.00295235387583076 -0.004887499696737629 … -0.0008444222562311587 0.002372371867952583; -0.008294524720988461 0.013731242514756659 … 0.002372371867952583 -0.006665087565281169]], Not converged after 5001 iterations.)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "include(\"testgaussianmle.jl\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.5.1", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.1" } }, "nbformat": 4, "nbformat_minor": 4 }