{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reduction to Hessenberg form\n",
    "\n",
    "The upper-triangular Hessenberg form has zeros below the first sub-diagonal\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "* & * & * & * & * & * \\\\\n",
    "* & * & * & * & * & * \\\\\n",
    "0 & * & * & * & * & * \\\\\n",
    "0 & 0 & * & * & * & * \\\\\n",
    "0 & 0 & 0 & * & * & * \\\\\n",
    "0 & 0 & 0 & 0 & * & * \\\\\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need sign function with sign(0) = 1\n",
    "def sign(x):\n",
    "    if x >= 0.0:\n",
    "        return 1.0\n",
    "    else:\n",
    "        return -1.0\n",
    "\n",
    "# Algorithm 26.1\n",
    "def hessenberg(A):\n",
    "    m = A.shape[0]\n",
    "    for k in range(m-2):\n",
    "        x = A[k+1:m, k]\n",
    "        e1 = np.zeros_like(x)\n",
    "        e1[0] = 1.0\n",
    "        v = sign(x[0]) * np.linalg.norm(x) * e1 + x\n",
    "        v = v / np.linalg.norm(v)\n",
    "        A[k+1:m,k:m] = A[k+1:m,k:m] - 2.0 * np.outer(v, np.dot(v, A[k+1:m,k:m]))\n",
    "        A[0:m,k+1:m] = A[0:m,k+1:m] - 2.0 * np.outer(A[0:m,k+1:m] @ v, v)\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test on a random matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.46e-02 -9.73e-01  3.66e-01 -2.42e-01  4.67e-01 -6.16e-01  6.15e-02]\n",
      " [-1.24e+00  2.15e+00 -8.66e-01  4.98e-03 -2.90e-01 -2.23e-01 -1.51e-01]\n",
      " [ 0.00e+00 -1.76e+00  1.74e-01 -1.22e-01  5.24e-01 -8.45e-03  1.60e-01]\n",
      " [-2.78e-17  2.22e-16  5.69e-01 -1.04e-02  1.19e-01  2.12e-01 -3.13e-01]\n",
      " [ 0.00e+00  2.22e-16  0.00e+00  2.86e-01  1.78e-01  1.56e-01 -5.46e-01]\n",
      " [ 0.00e+00 -2.22e-16  0.00e+00  0.00e+00  5.48e-02  1.30e-01 -1.65e-01]\n",
      " [ 0.00e+00 -8.67e-19  0.00e+00  0.00e+00 -8.67e-19  1.04e-01 -1.25e-01]]\n"
     ]
    }
   ],
   "source": [
    "m = 7\n",
    "A = np.random.rand(m,m)\n",
    "H = hessenberg(A)\n",
    "print(np.array_str(H, precision=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the lower triangular part which is supposed to be zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00]\n",
      " [ 0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00]\n",
      " [ 0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00]\n",
      " [-2.78e-17  2.22e-16  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00]\n",
      " [ 0.00e+00  2.22e-16  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00]\n",
      " [ 0.00e+00 -2.22e-16  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00]\n",
      " [ 0.00e+00 -8.67e-19  0.00e+00  0.00e+00 -8.67e-19  0.00e+00  0.00e+00]]\n"
     ]
    }
   ],
   "source": [
    "L = np.tril(H,k=-2)\n",
    "print(np.array_str(L, precision=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the maximum over the zero elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.220446049250313e-16\n"
     ]
    }
   ],
   "source": [
    "print(np.abs(L).max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test on a larger matrix without printing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1102230246251565e-15\n"
     ]
    }
   ],
   "source": [
    "m = 100\n",
    "A = np.random.rand(m, m)\n",
    "H = hessenberg(A)\n",
    "L = np.tril(H,k=-2)\n",
    "print(np.abs(L).max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hermitian case\n",
    "\n",
    "In this case, the Hessenberg reduction should give a tridiagonal matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.91 -1.34  0.    0.   -0.    0.   -0.  ]\n",
      " [-1.34  2.56 -1.57 -0.   -0.    0.   -0.  ]\n",
      " [ 0.   -1.57  0.79  0.31  0.   -0.   -0.  ]\n",
      " [ 0.   -0.    0.31 -0.09  0.26 -0.    0.  ]\n",
      " [ 0.    0.    0.    0.26 -0.21  0.23  0.  ]\n",
      " [ 0.    0.    0.    0.    0.23  0.41 -0.51]\n",
      " [-0.    0.    0.    0.   -0.   -0.51  0.03]]\n"
     ]
    }
   ],
   "source": [
    "m = 7\n",
    "B = np.random.rand(m,m)\n",
    "A = 0.5 * (B + B.T)\n",
    "H = hessenberg(A)\n",
    "print(np.array_str(H, precision=2, suppress_small=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remark\n",
    "\n",
    "The above function may not work correctly if the input matrix is complex. Write a complex version and test it."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  },
  "vscode": {
   "interpreter": {
    "hash": "c6e4e9f98eb68ad3b7c296f83d20e6de614cb42e90992a65aa266555a3137d0d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}