commit 54100fa836f12ad966578de13e8272fb02ae12b7
parent 1d9d6d8a7d80be6cd48a39b585318d5930eabdda
Author: miksa234 <milutin@popovic.xyz>
Date: Mon, 23 May 2022 07:12:12 +0200
foo:
Diffstat:
2 files changed, 205 insertions(+), 0 deletions(-)
diff --git a/num_ana/prog/.ipynb_checkpoints/prb9_p-checkpoint.ipynb b/num_ana/prog/.ipynb_checkpoints/prb9_p-checkpoint.ipynb
@@ -0,0 +1,6 @@
+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/num_ana/prog/prb9_p.ipynb b/num_ana/prog/prb9_p.ipynb
@@ -0,0 +1,199 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "0ae2f006",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 184,
+ "id": "6946cd07",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def fourier_frame(n: int, d: int) -> (np.ndarray, np.ndarray):\n",
+ " ω = lambda i, j : np.exp(-2*np.pi *1j/n * (i-1) * (j-1))\n",
+ " H = np.zeros([n, d], dtype='complex')\n",
+ " for i in range(n):\n",
+ " for j in range(d):\n",
+ " H[i, j] = 1/np.sqrt(n) * ω(i, j)\n",
+ " S = np.conjugate(H.T) @ H \n",
+ " return H, S \n",
+ "\n",
+ "def reconstruct(v: np.array, S: np.ndarray, H: np.ndarray) -> np.array:\n",
+ " s = 0\n",
+ " for j in range(H.shape[0]):\n",
+ " s += v @ np.conjugate(H.T)[:,j] * invs(S) @ H[j]\n",
+ " return np.round(s, 2)\n",
+ "\n",
+ "def qr_eig(A: np.ndarray, eps: float=1e-5) -> np.ndarray:\n",
+ " A_k = np.copy(A)\n",
+ " while True:\n",
+ " Q, R = np.linalg.qr(A_k)\n",
+ " A_1 = np.copy(A_k)\n",
+ " A_k = R@Q\n",
+ " if np.linalg.norm(A - A_1)/np.linalg.norm(A_1) < eps:\n",
+ " break\n",
+ " return np.diagonal(A_k)\n",
+ "\n",
+ "# in our particular case them matrix S is unitary\n",
+ "def invs(A: np.ndarray) -> np.ndarray:\n",
+ " return np.linalg.inv(A)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 185,
+ "id": "75be68d5",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "(1-0j) (1-0j) (1+0j) (1+0j) (1+0j) (1-0j)\n"
+ ]
+ }
+ ],
+ "source": [
+ "d = 6; n = 10\n",
+ "H, S = fourier_frame(n, d)\n",
+ "print(*np.round(qr_eig(S), 2))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 186,
+ "id": "fbe1004e",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "(1-0j) (2-0j) (3+0j) (4+0j) (5+0j) (6-0j)\n"
+ ]
+ }
+ ],
+ "source": [
+ "v = np.arange(1, d+1) \n",
+ "print(*reconstruct(v, S, H))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 187,
+ "id": "24f82555",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "(1.44-1.29j) (1.56-1.48j) (2.29+0.38j) (4+1.71j) (5.71+0.68j) (6.44-1.29j)\n"
+ ]
+ }
+ ],
+ "source": [
+ "H̃ = H[:-1]\n",
+ "S̃ = np.conjugate(H̃.T) @ H̃ \n",
+ "print(*reconstruct(v, S̃, H̃))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c5bdeb99",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "H, S = fourier_frame(5, 5)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "id": "ba0f89b5",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 1.+0.j, 0.-0.j, 0.+0.j, -0.+0.j, -0.-0.j],\n",
+ " [ 0.+0.j, 1.+0.j, 0.-0.j, 0.+0.j, -0.+0.j],\n",
+ " [ 0.-0.j, 0.+0.j, 1.+0.j, 0.-0.j, 0.+0.j],\n",
+ " [-0.-0.j, 0.-0.j, 0.+0.j, 1.+0.j, 0.-0.j],\n",
+ " [-0.+0.j, -0.-0.j, 0.-0.j, 0.+0.j, 1.+0.j]])"
+ ]
+ },
+ "execution_count": 14,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "np.round(np.conjugate(H.T) @ H, 2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "85380980",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "array([[ 0.-0.j, 0.+0.j, 1.+0.j, 0.-0.j, 0.+0.j],\n",
+ " [ 0.+0.j, 1.+0.j, 0.-0.j, 0.+0.j, -0.+0.j],\n",
+ " [ 1.+0.j, 0.-0.j, 0.+0.j, -0.+0.j, -0.-0.j],\n",
+ " [ 0.-0.j, 0.+0.j, -0.+0.j, -0.-0.j, 1.+0.j],\n",
+ " [ 0.+0.j, -0.+0.j, -0.-0.j, 1.+0.j, 0.-0.j]])"
+ ]
+ },
+ "execution_count": 19,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "np.round(H@H, 2)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4e2bedda",
+ "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.10.4"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
