hw1_conv1D.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "053e3f07",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4587586843129764\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "class Convolution1d :\n",
    "    def __init__(self, filt) :\n",
    "        self.__filt = filt\n",
    "        self.__r = filt.size\n",
    "        self.T = TransposedConvolution1d(self.__filt)\n",
    "\n",
    "    def __matmul__(self, vector) :\n",
    "        r, n = self.__r, vector.size\n",
    "        \n",
    "        return np.asarray([ np.dot(self.__filt, vector[i:i+r]) for i in range(n-r+1) ])\n",
    "\n",
    "\n",
    "class TransposedConvolution1d :\n",
    "    \"\"\"\n",
    "    Transpose of 1-dimensional convolution operator used for the\n",
    "    transpose-convolution operation A.T@(...)\n",
    "    \"\"\"\n",
    "    def __init__(self, filt) :\n",
    "        self.__filt = filt\n",
    "        self.__r = filt.size\n",
    "\n",
    "    def __matmul__(self, vector) :\n",
    "        r = self.__r\n",
    "        n = vector.size + r - 1\n",
    "\n",
    "        vector_prime = np.concatenate((np.zeros(r-1), vector, np.zeros(r-1)))\n",
    "        \n",
    "        return np.asarray( [np.dot(self.__filt[::-1], vector_prime[i:i+r]) for i in range(n)] )\n",
    "\n",
    "def huber_loss(x) :\n",
    "    return np.sum( (1/2)*(x**2)*(np.abs(x)<=1) + (np.sign(x)*x-1/2)*(np.abs(x)>1) )\n",
    "\n",
    "def huber_grad(x) :\n",
    "    return x*(np.abs(x)<=1) + np.sign(x)*(np.abs(x)>1)\n",
    "\n",
    "r, n, lam = 3, 20, 0.1\n",
    "\n",
    "np.random.seed(0)\n",
    "k = np.random.randn(r)\n",
    "b = np.random.randn(n-r+1)\n",
    "\n",
    "A = Convolution1d(k)\n",
    "#from scipy.linalg import circulant\n",
    "#A = circulant(np.concatenate((np.flip(k),np.zeros(n-r))))[r-1:,:]\n",
    "\n",
    "x = np.zeros(n)\n",
    "alpha = 0.01\n",
    "for _ in range(100) :\n",
    "    x = x - alpha*(A.T@(huber_grad(A@x-b))+lam*x)\n",
    "\n",
    "print(huber_loss(A@x-b)+0.5*lam*np.linalg.norm(x)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "375fb840",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}