Test script of the positive semidefinite projection of a symmetric matrix.

numerical_error.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This is the Julia code that I use for the svd.\n",
    "\n",
    "**Create a random symmetric matrix:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 0.590845  0.854147   0.648882   0.112486   0.950498\n",
       " 0.854147  0.200586   0.0109059  0.276021   0.96467 \n",
       " 0.648882  0.0109059  0.066423   0.651664   0.945775\n",
       " 0.112486  0.276021   0.651664   0.0566425  0.789904\n",
       " 0.950498  0.96467    0.945775   0.789904   0.82116 "
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Base.Test\n",
    "rng = MersenneTwister(1234);\n",
    "X = rand(rng,5,5);\n",
    "X = full(Symmetric(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compute eigenvalue decomposition:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues: [-1.01743, -0.654582, -0.108438, 0.445086, 3.07102]\n",
      "Eigenvectors:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       "  0.331735  -0.35125     0.554806  -0.47478   -0.483062\n",
       " -0.553521  -0.167774   -0.528168  -0.488244  -0.384867\n",
       " -0.638086   0.0619594   0.493018   0.458838  -0.367978\n",
       "  0.242374  -0.653637   -0.342732   0.555047  -0.297439\n",
       "  0.343026   0.646061   -0.229555   0.132634  -0.628213"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = eigfact(X)\n",
    "Λ = diagm(F[:values])\n",
    "Q = F[:vectors]\n",
    "println(\"Eigenvalues: $(F[:values])\")\n",
    "println(\"Eigenvectors:\")\n",
    "Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "X has 3 negative eigenvalues.\n",
    "\n",
    "**Perform projection onto sd cone:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 0.816949  0.674123  0.448932  0.323957  0.903922\n",
       " 0.674123  0.560988  0.335215  0.230936  0.713684\n",
       " 0.448932  0.335215  0.509545  0.44948   0.73701 \n",
       " 0.323957  0.230936  0.44948   0.408814  0.606602\n",
       " 0.903922  0.713684  0.73701   0.606602  1.21981 "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xp = Q*max.(Λ,0)*Q'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Check the eigenvalues of the projected matrix Xp:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues: [3.07102, -2.63735e-16, 0.445086, -2.65217e-17, 7.12942e-17]\n"
     ]
    }
   ],
   "source": [
    "E = eigfact(Xp)\n",
    "println(\"Eigenvalues: $(E[:values])\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Through numerical errors, Xp has some negative eigenvalues very close to zero.\n",
    "\n",
    "**Take an eigenvector v corresponding to an almost zero eigenvalue, i.e. the second one**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = E[:vectors][:,2];\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Compare the following according to Paul:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.708201408659274e-17"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = v'*Xp*v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1606330598889774e-32"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = (v'*Q)*max.(Λ,0)*(Q'*v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[91mTest Failed\n",
      "\u001b[39m\u001b[22m  Expression: A == B\n",
      "   Evaluated: -1.708201408659274e-17 == 0.2604210834761241\n"
     ]
    },
    {
     "ename": "LoadError",
     "evalue": "\u001b[91mThere was an error during testing\u001b[39m",
     "output_type": "error",
     "traceback": [
      "\u001b[91mThere was an error during testing\u001b[39m",
      "",
      "Stacktrace:",
      " [1] \u001b[1mrecord\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Base.Test.FallbackTestSet, ::Base.Test.Fail\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./test.jl:533\u001b[22m\u001b[22m",
      " [2] \u001b[1mdo_test\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::Base.Test.Returned, ::Expr\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./test.jl:352\u001b[22m\u001b[22m",
      " [3] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:515\u001b[22m\u001b[22m"
     ]
    }
   ],
   "source": [
    "@test A == B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two expressions are obviously not the same, so it seems like the error is introduced by the second refactorization!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}