Let $A$ be the $20000 \times 20000$ matrix whose entries are 0 everywhere except for the primes 2,3,5,7,...,224737 along the main diagonal and the number 1 in all the positions $ a_{ij} $ with $ |i-j| = 1,2,4,8,\cdots,16384$. What is the (1,1) entry of $A^{-1}$
def is_prime(N):
if N < 2: return False
if N == 2: return True
for i in range(2, int(N**0.5)+1):
if N%i == 0: return False
return True
def is_prime(N):
return N == 2 or any([N%i == 0 for i in range()])
entries = []
primes = [i for i in range(2,250000) if is_prime(i)][:20000]
entries.extend([(i,i,primes[i]) for i in range(len(primes))])
for diff in [2**j for j in range(15)]:
for i in range(20000):
j1, j2 = i+diff, i-diff
if 0<=j1<20000: entries.append((i,j1,1))
if 0<=j2<20000: entries.append((i,j2,1))
print(len(entries))
554466
We avoid calculating $A^{-1}$ explicitly. Instead, note that if we have $\hat{v}$ be the first column of $ A^{-1} $, then we get that $ A \hat{v} = (1,0,\cdots, 0)^{T}$. Let's solve this system.
import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve
row = np.array([i[0] for i in entries])
col = np.array([i[1] for i in entries])
data = np.array([i[2] for i in entries])
A = coo_matrix((data, (row, col))).tocsc()
b = np.zeros(A.shape[0])
b[0] = 1
v = spsolve(A,b)
# Timeout
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-e0aa61c5c910> in <module>() 3 from scipy.sparse.linalg import spsolve 4 ----> 5 row = np.array([i[0] for i in entries]) 6 col = np.array([i[1] for i in entries]) 7 data = np.array([i[2] for i in entries]) NameError: name 'entries' is not defined
The above code times out - so instead, let's try to use an iterative method to solve this linear system.
import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import bicg
row = np.array([i[0] for i in entries])
col = np.array([i[1] for i in entries])
data = np.array([i[2] for i in entries])
A = coo_matrix((data, (row, col))).tocsc()
b = np.zeros(A.shape[0])
b[0] = 1
iterations = 0
def callback(xk):
global iterations
iterations += 1
if iterations % 10 == 0: print("Iteration {}: current val = {}".format(iterations, xk[0]))
x, exit_code = bicg(A,b,tol=1e-50,callback=callback)
print(x[0])
Iteration 10: current val = 0.5160474805069852 Iteration 20: current val = 0.5424188205227113 Iteration 30: current val = 0.5649984704641765 Iteration 40: current val = 0.5809367071267324 Iteration 50: current val = 0.6096730237142656 Iteration 60: current val = 0.6190393486130085 Iteration 70: current val = 0.6263245512238234 Iteration 80: current val = 0.6508621776392641 Iteration 90: current val = 0.6577928020000656 Iteration 100: current val = 0.6608853920165397 Iteration 110: current val = 0.6646370853741161 Iteration 120: current val = 0.6769767188475972 Iteration 130: current val = 0.686598661153931 Iteration 140: current val = 0.6950283972972651 Iteration 150: current val = 0.6985776576929462 Iteration 160: current val = 0.700549226696957 Iteration 170: current val = 0.7016993660417654 Iteration 180: current val = 0.7037530670049756 Iteration 190: current val = 0.70578202448334 Iteration 200: current val = 0.7080325880243894 Iteration 210: current val = 0.7107390862052434 Iteration 220: current val = 0.7131266080608784 Iteration 230: current val = 0.7154864307143483 Iteration 240: current val = 0.7169760206590395 Iteration 250: current val = 0.7177222840589774 Iteration 260: current val = 0.7182998666211887 Iteration 270: current val = 0.7189625336279692 Iteration 280: current val = 0.7195270338881113 Iteration 290: current val = 0.7198970070328818 Iteration 300: current val = 0.7202715787882576 Iteration 310: current val = 0.720548747555294 Iteration 320: current val = 0.7208623428688858 Iteration 330: current val = 0.7211885297932319 Iteration 340: current val = 0.7216667311611941 Iteration 350: current val = 0.7219940973832849 Iteration 360: current val = 0.7224389483334973 Iteration 370: current val = 0.7228169933348944 Iteration 380: current val = 0.7231983252626363 Iteration 390: current val = 0.7234404924750876 Iteration 400: current val = 0.7236478043297795 Iteration 410: current val = 0.7238808716336398 Iteration 420: current val = 0.7241166500666772 Iteration 430: current val = 0.7243659045725263 Iteration 440: current val = 0.7244940472057783 Iteration 450: current val = 0.724633608525676 Iteration 460: current val = 0.7247020541105452 Iteration 470: current val = 0.7247680447799154 Iteration 480: current val = 0.7248237539591655 Iteration 490: current val = 0.7248677969693484 Iteration 500: current val = 0.7248916848307266 Iteration 510: current val = 0.7249208479337073 Iteration 520: current val = 0.7249533593512273 Iteration 530: current val = 0.7249736254056114 Iteration 540: current val = 0.7249887123438771 Iteration 550: current val = 0.7250067816490801 Iteration 560: current val = 0.7250207131481929 Iteration 570: current val = 0.7250304954206683 Iteration 580: current val = 0.7250427286641773 Iteration 590: current val = 0.72505125564296 Iteration 600: current val = 0.7250561238493789 Iteration 610: current val = 0.7250605939852699 Iteration 620: current val = 0.7250640020707945 Iteration 630: current val = 0.7250660938618454 Iteration 640: current val = 0.7250677953689088 Iteration 650: current val = 0.7250693490905185 Iteration 660: current val = 0.725070446117053 Iteration 670: current val = 0.7250713225273709 Iteration 680: current val = 0.7250726985038631 Iteration 690: current val = 0.7250739362881828 Iteration 700: current val = 0.7250748779230066 Iteration 710: current val = 0.7250755450783181 Iteration 720: current val = 0.7250760283885833 Iteration 730: current val = 0.7250765020947028 Iteration 740: current val = 0.7250770006680382 Iteration 750: current val = 0.7250773496606184 Iteration 760: current val = 0.7250776006228311 Iteration 770: current val = 0.7250777636744249 Iteration 780: current val = 0.7250778756720382 Iteration 790: current val = 0.7250779688442739 Iteration 800: current val = 0.7250780439687791 Iteration 810: current val = 0.7250781067493376 Iteration 820: current val = 0.7250781516265811 Iteration 830: current val = 0.7250781883759363 Iteration 840: current val = 0.7250782241898078 Iteration 850: current val = 0.7250782614326434 Iteration 860: current val = 0.7250782837700643 Iteration 870: current val = 0.7250782994297919 Iteration 880: current val = 0.7250783112476437 Iteration 890: current val = 0.7250783221719436 Iteration 900: current val = 0.725078329800566 Iteration 910: current val = 0.7250783351385289 Iteration 920: current val = 0.7250783385028293 Iteration 930: current val = 0.7250783410043907 Iteration 940: current val = 0.725078342638123 Iteration 950: current val = 0.7250783438331366 Iteration 960: current val = 0.7250783446050801 Iteration 970: current val = 0.7250783450804342 Iteration 980: current val = 0.7250783453953678 Iteration 990: current val = 0.725078345614237 Iteration 1000: current val = 0.7250783457924285 Iteration 1010: current val = 0.7250783459042496 Iteration 1020: current val = 0.7250783459839163 Iteration 1030: current val = 0.7250783460699098 Iteration 1040: current val = 0.7250783461385248 Iteration 1050: current val = 0.7250783461804372 Iteration 1060: current val = 0.7250783462047834 Iteration 1070: current val = 0.7250783462184939 Iteration 1080: current val = 0.7250783462318636 Iteration 1090: current val = 0.7250783462415041 Iteration 1100: current val = 0.7250783462487344 Iteration 1110: current val = 0.7250783462541678 Iteration 1120: current val = 0.7250783462581615 Iteration 1130: current val = 0.72507834626106 Iteration 1140: current val = 0.725078346263089 Iteration 1150: current val = 0.72507834626454 Iteration 1160: current val = 0.7250783462657099 Iteration 1170: current val = 0.7250783462665752 Iteration 1180: current val = 0.7250783462671782 Iteration 1190: current val = 0.7250783462676181 Iteration 1200: current val = 0.7250783462679155 Iteration 1210: current val = 0.7250783462681067 Iteration 1220: current val = 0.7250783462682162 Iteration 1230: current val = 0.7250783462682798 Iteration 1240: current val = 0.7250783462683179 Iteration 1250: current val = 0.725078346268345 Iteration 1260: current val = 0.7250783462683673 Iteration 1270: current val = 0.7250783462683824 Iteration 1280: current val = 0.7250783462683897 Iteration 1290: current val = 0.7250783462683944 Iteration 1300: current val = 0.7250783462683967 Iteration 1310: current val = 0.7250783462683981 Iteration 1320: current val = 0.7250783462683992 Iteration 1330: current val = 0.7250783462684001 Iteration 1340: current val = 0.7250783462684001 Iteration 1350: current val = 0.7250783462684001 Iteration 1360: current val = 0.7250783462684001 Iteration 1370: current val = 0.7250783462684001 Iteration 1380: current val = 0.7250783462684001 Iteration 1390: current val = 0.7250783462684001 Iteration 1400: current val = 0.7250783462684001 Iteration 1410: current val = 0.7250783462684001 Iteration 1420: current val = 0.7250783462684001 Iteration 1430: current val = 0.7250783462684001 Iteration 1440: current val = 0.7250783462684001 Iteration 1450: current val = 0.7250783462684001 Iteration 1460: current val = 0.7250783462684001 Iteration 1470: current val = 0.7250783462684001 Iteration 1480: current val = 0.7250783462684001 Iteration 1490: current val = 0.7250783462684001 Iteration 1500: current val = 0.7250783462684001 Iteration 1510: current val = 0.7250783462684001 Iteration 1520: current val = 0.7250783462684001 Iteration 1530: current val = 0.7250783462684001 Iteration 1540: current val = 0.7250783462684001 Iteration 1550: current val = 0.7250783462684001 Iteration 1560: current val = 0.7250783462684001 Iteration 1570: current val = 0.7250783462684001 Iteration 1580: current val = 0.7250783462684001 Iteration 1590: current val = 0.7250783462684001 Iteration 1600: current val = 0.7250783462684001 Iteration 1610: current val = 0.7250783462684001 Iteration 1620: current val = 0.7250783462684001 Iteration 1630: current val = 0.7250783462684001 Iteration 1640: current val = 0.7250783462684001 Iteration 1650: current val = 0.7250783462684001 Iteration 1660: current val = 0.7250783462684001 Iteration 1670: current val = 0.7250783462684001 Iteration 1680: current val = 0.7250783462684001 Iteration 1690: current val = 0.7250783462684001 Iteration 1700: current val = 0.7250783462684001 Iteration 1710: current val = 0.7250783462684001 Iteration 1720: current val = 0.7250783462684001 Iteration 1730: current val = 0.7250783462684001 Iteration 1740: current val = 0.7250783462684001 Iteration 1750: current val = 0.7250783462684001 Iteration 1760: current val = 0.7250783462684001 Iteration 1770: current val = 0.7250783462684001 Iteration 1780: current val = 0.7250783462684001 Iteration 1790: current val = 0.7250783462684001 Iteration 1800: current val = 0.7250783462684001 Iteration 1810: current val = 0.7250783462684001 Iteration 1820: current val = 0.7250783462684001 Iteration 1830: current val = 0.7250783462684001 Iteration 1840: current val = 0.7250783462684001 Iteration 1850: current val = 0.7250783462684001 Iteration 1860: current val = 0.7250783462684001 Iteration 1870: current val = 0.7250783462684001 Iteration 1880: current val = 0.7250783462684001 Iteration 1890: current val = 0.7250783462684001 Iteration 1900: current val = 0.7250783462684001 Iteration 1910: current val = 0.7250783462684001 Iteration 1920: current val = 0.7250783462684001 Iteration 1930: current val = 0.7250783462684001 Iteration 1940: current val = 0.7250783462684001 Iteration 1950: current val = 0.7250783462684001 Iteration 1960: current val = 0.7250783462684001 Iteration 1970: current val = 0.7250783462684001 Iteration 1980: current val = 0.7250783462684001 Iteration 1990: current val = 0.7250783462684001 Iteration 2000: current val = 0.7250783462684001 Iteration 2010: current val = 0.7250783462684001 Iteration 2020: current val = 0.7250783462684001 Iteration 2030: current val = 0.7250783462684001 Iteration 2040: current val = 0.7250783462684001 Iteration 2050: current val = 0.7250783462684001 Iteration 2060: current val = 0.7250783462684001 Iteration 2070: current val = 0.7250783462684001 Iteration 2080: current val = 0.7250783462684001 Iteration 2090: current val = 0.7250783462684001 Iteration 2100: current val = 0.7250783462684001 Iteration 2110: current val = 0.7250783462684001 Iteration 2120: current val = 0.7250783462684001 Iteration 2130: current val = 0.7250783462684001 Iteration 2140: current val = 0.7250783462684001 Iteration 2150: current val = 0.7250783462684001 Iteration 2160: current val = 0.7250783462684001 Iteration 2170: current val = 0.7250783462684001 Iteration 2180: current val = 0.7250783462684001 Iteration 2190: current val = 0.7250783462684001 Iteration 2200: current val = 0.7250783462684001 Iteration 2210: current val = 0.7250783462684001 Iteration 2220: current val = 0.7250783462684001 Iteration 2230: current val = 0.7250783462684001 Iteration 2240: current val = 0.7250783462684001 Iteration 2250: current val = 0.7250783462684001 Iteration 2260: current val = 0.7250783462684001 Iteration 2270: current val = 0.7250783462684001 Iteration 2280: current val = 0.7250783462684001 Iteration 2290: current val = 0.7250783462684001 0.7250783462684001