Skip to content

Commit d7063c6

Browse files
lorenzotomadalorenzotomada
committed
updating QR
1 parent 542cf89 commit d7063c6

File tree

4 files changed

+251
-125
lines changed

4 files changed

+251
-125
lines changed

‎docs/Documentation.ipynb

Lines changed: 84 additions & 74 deletions
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,7 @@
7373
"from line_profiler import profile, LineProfiler\n",
7474
"from time import time\n",
7575
"\n",
76+
"\n",
7677
"def Lanczos_PRO_original(A, q, m=None, toll=np.sqrt(np.finfo(float).eps)):\n",
7778
" \"\"\"\n",
7879
" Perform the Lanczos algorithm for symmetric matrices.\n",
@@ -134,12 +135,11 @@
134135
" if np.abs(beta[j]) < 1e-15:\n",
135136
" return Q, alpha, beta[:-1]\n",
136137
"\n",
137-
"\n",
138138
" return Q, alpha, beta[:-1]\n",
139139
"\n",
140140
"\n",
141-
"#@jit(parallel=True) \n",
142-
"def Lanczos_PRO(A, q, m=None, toll=np.sqrt(np.finfo(float).eps)):\n",
141+
"# @jit(parallel=True)\n",
142+
"def Lanczos_PRO(A, q, m=None, toll=np.sqrt(np.finfo(float).eps)):\n",
143143
" \"\"\"\n",
144144
" Lanczos algorithm for symmetric matrices\n",
145145
" :param A: symmetric matrix square matrix of size n\n",
@@ -152,40 +152,39 @@
152152
" alpha: vector of size m. Diagonal of the tridiagonal matrix\n",
153153
" beta: vector of size m-1. Upper diagonal of the tridiagonal matrix\n",
154154
" \"\"\"\n",
155-
" if m==None:\n",
156-
" m=A.shape[0]\n",
157-
" \n",
158-
" q=q/np.linalg.norm(q)\n",
159-
" #Q=np.array([q])\n",
155+
" if m == None:\n",
156+
" m = A.shape[0]\n",
157+
"\n",
158+
" q = q / np.linalg.norm(q)\n",
159+
" # Q=np.array([q])\n",
160160
" Q = np.zeros((m, A.shape[0]))\n",
161161
" Q[0] = q\n",
162-
" r=A@q\n",
163-
" alpha=[]\n",
164-
" beta=[]\n",
165-
" alpha.append(q@r)\n",
166-
" r=r-alpha[0]*q\n",
162+
" r = A @ q\n",
163+
" alpha = []\n",
164+
" beta = []\n",
165+
" alpha.append(q @ r)\n",
166+
" r = r - alpha[0] * q\n",
167167
" beta.append(np.linalg.norm(r))\n",
168168
"\n",
169169
" for j in range(1, m):\n",
170-
" q=r/beta[j-1]\n",
171-
" if np.any(np.abs(q@Q[:j-1].T)>toll):\n",
172-
" for q_bbasis in Q[:j-1]:\n",
170+
" q = r / beta[j - 1]\n",
171+
" if np.any(np.abs(q @ Q[: j - 1].T) > toll):\n",
172+
" for q_bbasis in Q[: j - 1]:\n",
173173
" q = q - (q @ q_bbasis) * q_bbasis\n",
174174
" # for i in prange(j-1):\n",
175175
" # q = q - (q @ Q[i]) * Q[i]\n",
176176
"\n",
177-
" q=q/np.linalg.norm(q)\n",
178-
" Q[j]=q\n",
179-
" r=A@q-beta[j-1]*Q[j-1]\n",
180-
" alpha.append(q@r)\n",
181-
" r=r-alpha[j]*q\n",
177+
" q = q / np.linalg.norm(q)\n",
178+
" Q[j] = q\n",
179+
" r = A @ q - beta[j - 1] * Q[j - 1]\n",
180+
" alpha.append(q @ r)\n",
181+
" r = r - alpha[j] * q\n",
182182
" beta.append(np.linalg.norm(r))\n",
183183
"\n",
184-
" if np.abs(beta[j])<1e-15:\n",
185-
" \n",
184+
" if np.abs(beta[j]) < 1e-15:\n",
185+
"\n",
186186
" return Q, alpha, beta[:-1]\n",
187-
" return Q, alpha, beta[:-1]\n",
188-
"\n"
187+
" return Q, alpha, beta[:-1]"
189188
]
190189
},
191190
{
@@ -331,24 +330,22 @@
331330
"size = 1000\n",
332331
"A = np.random.rand(size, size)\n",
333332
"A = (A + A.T) / 2\n",
334-
"A=np.array(A, dtype=float)\n",
333+
"A = np.array(A, dtype=float)\n",
335334
"q = np.random.rand(size)\n",
336335
"lp = LineProfiler()\n",
337336
"lp.add_function(Lanczos_PRO_original)\n",
338-
"result = lp.run('Lanczos_PRO_original(A, q)')\n",
339-
"#lp.print_stats()\n",
340-
"#Q, alpha, beta = Lanczos_PRO(A, q)\n",
337+
"result = lp.run(\"Lanczos_PRO_original(A, q)\")\n",
338+
"# lp.print_stats()\n",
339+
"# Q, alpha, beta = Lanczos_PRO(A, q)\n",
341340
"\n",
342341
"\n",
343-
"\n",
344-
"t_s=time()\n",
342+
"t_s = time()\n",
345343
"lp.add_function(Lanczos_PRO)\n",
346-
"result = lp.run('Lanczos_PRO(A, q)')\n",
344+
"result = lp.run(\"Lanczos_PRO(A, q)\")\n",
347345
"lp.print_stats()\n",
348346
"\n",
349-
"t_e=time()\n",
350-
"print(f\"Optimized function time: {t_e-t_s}\")\n",
351-
"\n"
347+
"t_e = time()\n",
348+
"print(f\"Optimized function time: {t_e-t_s}\")"
352349
]
353350
},
354351
{
@@ -388,7 +385,7 @@
388385
" A = A_copy.copy()\n",
389386
" iter = 0\n",
390387
" Q = np.eye(A.shape[0])\n",
391-
" \n",
388+
"\n",
392389
" # Correctly preallocate as a 2D array (n-1, 2)\n",
393390
" Matrix_trigonometry = np.zeros((A.shape[0] - 1, 2))\n",
394391
"\n",
@@ -397,22 +394,31 @@
397394
" for i in range(A.shape[0] - 1):\n",
398395
" c = A[i, i] / np.sqrt(A[i, i] ** 2 + A[i + 1, i] ** 2)\n",
399396
" s = -A[i + 1, i] / np.sqrt(A[i, i] ** 2 + A[i + 1, i] ** 2)\n",
400-
" Matrix_trigonometry[i, :] = [c, s] \n",
401-
" R=np.eye(A.shape[0])\n",
397+
" Matrix_trigonometry[i, :] = [c, s]\n",
398+
" R = np.eye(A.shape[0])\n",
402399
" # Apply the Givens rotation to A (modify in place)\n",
403-
" R[i:i+2, i:i+2] = np.array([[c, -s], [s, c]])\n",
404-
" A= R @ A\n",
405-
" A[i+1, i] = 0 \n",
406-
"\n",
407-
"\n",
408-
" Q=np.eye(A.shape[0])\n",
409-
" i=0\n",
410-
" Q[0:2, 0:2]=np.array([[Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]], [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]]])\n",
411-
" for i in range(1, A.shape[0]-1):\n",
412-
" R=np.eye(A.shape[0])\n",
413-
" R[i: i+2, i:i+2]=np.array([[Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]], [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]]])\n",
414-
" Q = Q@R\n",
415-
" A=A@Q\n",
400+
" R[i : i + 2, i : i + 2] = np.array([[c, -s], [s, c]])\n",
401+
" A = R @ A\n",
402+
" A[i + 1, i] = 0\n",
403+
"\n",
404+
" Q = np.eye(A.shape[0])\n",
405+
" i = 0\n",
406+
" Q[0:2, 0:2] = np.array(\n",
407+
" [\n",
408+
" [Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]],\n",
409+
" [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]],\n",
410+
" ]\n",
411+
" )\n",
412+
" for i in range(1, A.shape[0] - 1):\n",
413+
" R = np.eye(A.shape[0])\n",
414+
" R[i : i + 2, i : i + 2] = np.array(\n",
415+
" [\n",
416+
" [Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]],\n",
417+
" [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]],\n",
418+
" ]\n",
419+
" )\n",
420+
" Q = Q @ R\n",
421+
" A = A @ Q\n",
416422
" A = A @ Q # Update A\n",
417423
" iter += 1\n",
418424
"\n",
@@ -421,43 +427,49 @@
421427
"\n",
422428
" return np.diag(A), Q\n",
423429
"\n",
424-
"#@jit(nopython=True)\n",
430+
"\n",
431+
"# @jit(nopython=True)\n",
425432
"def QR_method_optimized(A_copy, tol=1e-10, max_iter=100):\n",
426433
" A = A_copy.copy()\n",
427434
" iter = 0\n",
428435
" Q = np.eye(A.shape[0])\n",
429-
" \n",
436+
"\n",
430437
" # Correctly preallocate as a 2D array (n-1, 2)\n",
431438
" Matrix_trigonometry = np.zeros((A.shape[0] - 1, 2))\n",
432-
" d=np.zeros(A.shape[0])\n",
433-
" #while np.linalg.norm((np.diag(A, -1)), np.inf) > tol and iter < max_iter:\n",
434-
" while np.linalg.norm((np.diag(A, 0)-d)/np.diag(A, 0), np.inf) > tol and iter < max_iter:\n",
439+
" d = np.zeros(A.shape[0])\n",
440+
" # while np.linalg.norm((np.diag(A, -1)), np.inf) > tol and iter < max_iter:\n",
441+
" while (\n",
442+
" np.linalg.norm((np.diag(A, 0) - d) / np.diag(A, 0), np.inf) > tol\n",
443+
" and iter < max_iter\n",
444+
" ):\n",
435445
" # Compute Givens rotation\n",
436-
" d=np.diag(A, 0)\n",
446+
" d = np.diag(A, 0)\n",
437447
" for i in range(A.shape[0] - 1):\n",
438448
" c = A[i, i] / np.sqrt(A[i, i] ** 2 + A[i + 1, i] ** 2)\n",
439449
" s = -A[i + 1, i] / np.sqrt(A[i, i] ** 2 + A[i + 1, i] ** 2)\n",
440-
" Matrix_trigonometry[i, :] = [c, s] \n",
450+
" Matrix_trigonometry[i, :] = [c, s]\n",
441451
"\n",
442452
" # Apply the Givens rotation to A (modify in place)\n",
443453
" R = np.array([[c, -s], [s, c]])\n",
444-
" A[i:i+2, i:] = R @ A[i:i+2, i:]\n",
445-
" A[i+1, i] = 0 \n",
454+
" A[i : i + 2, i:] = R @ A[i : i + 2, i:]\n",
455+
" A[i + 1, i] = 0\n",
446456
"\n",
447457
" # Construct full Q matrix from stored Givens rotations\n",
448458
" Q = np.eye(A.shape[0])\n",
449459
" for i in range(A.shape[0] - 1):\n",
450-
" R=np.array([[Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]], [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]]])\n",
451-
" Q[:, i:i+2] = Q[:, i:i+2]@R\n",
460+
" R = np.array(\n",
461+
" [\n",
462+
" [Matrix_trigonometry[i, 0], Matrix_trigonometry[i, 1]],\n",
463+
" [-Matrix_trigonometry[i, 1], Matrix_trigonometry[i, 0]],\n",
464+
" ]\n",
465+
" )\n",
466+
" Q[:, i : i + 2] = Q[:, i : i + 2] @ R\n",
452467
" A = A @ Q # Update A\n",
453468
" iter += 1\n",
454469
"\n",
455470
" return np.diag(A), Q\n",
456471
"\n",
457472
"\n",
458-
"\n",
459-
"\n",
460-
"\n",
461473
"np.random.seed(0)\n",
462474
"size = 300\n",
463475
"A = np.random.rand(size, size)\n",
@@ -466,21 +478,20 @@
466478
"\n",
467479
"q, alpha, beta = Lanczos_PRO(A, x0, size)\n",
468480
"\n",
469-
"T=np.diag(alpha)+np.diag(beta, k=1)+np.diag(beta, k=-1)\n",
481+
"T = np.diag(alpha) + np.diag(beta, k=1) + np.diag(beta, k=-1)\n",
470482
"QR_method_optimized(T)\n",
471483
"lp_QR = LineProfiler()\n",
472484
"lp_QR.add_function(QR_method)\n",
473-
"result = lp_QR.run('QR_method(T)')\n",
485+
"result = lp_QR.run(\"QR_method(T)\")\n",
474486
"\n",
475487
"\n",
476488
"lp_QR.add_function(QR_method_optimized)\n",
477-
"result = lp_QR.run('QR_method_optimized(T)')\n",
489+
"result = lp_QR.run(\"QR_method_optimized(T)\")\n",
478490
"lp_QR.print_stats()\n",
479491
"# t_s=time()\n",
480492
"# result= QR_method_optimized(T)\n",
481493
"# t_e=time()\n",
482-
"# print(f\"Optimized function time: {t_e-t_s}\")\n",
483-
"\n"
494+
"# print(f\"Optimized function time: {t_e-t_s}\")"
484495
]
485496
},
486497
{
@@ -513,18 +524,17 @@
513524
"x0 = np.random.rand(size)\n",
514525
"q, alpha, beta = Lanczos_PRO(A, x0, size)\n",
515526
"\n",
516-
"t_s=time()\n",
527+
"t_s = time()\n",
517528
"np.linalg.eig(A)\n",
518-
"t_e=time()\n",
529+
"t_e = time()\n",
519530
"\n",
520531
"print(f\"{t_e-t_s: .4e}\")\n",
521532
"\n",
522533
"# T=np.diag(alpha)+np.diag(beta, k=1)+np.diag(beta, k=-1)\n",
523534
"# t_s=time()\n",
524535
"# result= QR_method_optimized(T, max_iter=5000, tol=1e-6)\n",
525536
"# t_e=time()\n",
526-
"# print(f\"Optimized function time: {t_e-t_s}\")\n",
527-
"\n"
537+
"# print(f\"Optimized function time: {t_e-t_s}\")"
528538
]
529539
}
530540
],

‎src/pyclassify/__init__.py

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
"power_method_cp",
88
"Lanczos_PRO",
99
"QR_method",
10+
"QR",
1011
]
1112

1213
from .eigenvalues import (
@@ -18,4 +19,5 @@
1819
power_method_cp,
1920
Lanczos_PRO,
2021
QR_method,
22+
QR,
2123
)

0 commit comments

Comments
 (0)