In [2]:
import numpy as np
from scipy.linalg import svd
from scipy.sparse.linalg import svds
import pprint
import math
#=====================
# define matrix
#=====================
A = np.array([[1.,2.,3.],
[4.,5.,6.],
[7.,8.,9.],
[10.,11.,12.]])
print(A)
U, S, V = svd(A, lapack_driver="gesvd")
[[ 1. 2. 3.] [ 4. 5. 6.] [ 7. 8. 9.] [10. 11. 12.]]
In [3]:
print("U matrix:", U)
print("V matrix:", V)
print("singular values :", S)
U matrix: [[-0.14087668 0.82471435 0.54212637 -0.07809606] [-0.34394629 0.42626394 -0.6646258 0.50820522] [-0.54701591 0.02781353 -0.29712753 -0.78212226] [-0.75008553 -0.37063688 0.41962695 0.3520131 ]] V matrix: [[-0.50453315 -0.5745157 -0.64449826] [-0.76077568 -0.05714052 0.64649464] [-0.40824829 0.81649658 -0.40824829]] singular values : [2.54624074e+01 1.29066168e+00 1.42996198e-15]
In [3]:
# print(S[0])
# print((U[:,0][np.newaxis]).T)
# print(V[0,:][np.newaxis])
In [4]:
A1 = S[0] * (U[:,0][np.newaxis]).T @ V[0,:][np.newaxis]
print(A1)
[[ 1.80979033 2.06082192 2.31185351] [ 4.41855027 5.03143657 5.64432287] [ 7.02731022 8.00205122 8.97679223] [ 9.63607016 10.97266587 12.30926159]]
In [5]:
A1Diff = A1 - A
print(A1Diff)
[[ 0.80979033 0.06082192 -0.68814649] [ 0.41855027 0.03143657 -0.35567713] [ 0.02731022 0.00205122 -0.02320777] [-0.36392984 -0.02733413 0.30926159]]
In [6]:
F_Norm = 0
for i in range(len(A1Diff)): #0:4
for j in range(len(A1Diff[0])): #0:3
F_Norm += A1Diff[i][j]**2
F_Norm = pow(F_Norm,.5)
print("F_Norm: ",F_Norm)
F_Norm: 1.2906616757612315
In [7]:
SingErr = 0
for i in range(1,len(S)):
SingErr += pow(S[i],2)
SingErr = pow(SingErr,.5)
print("Singular Err: ",SingErr)
Singular Err: 1.290661675761233
In [10]:
#2.
import matplotlib.pyplot as plt
import cv2
im = cv2.imread('boat.png')
imMatrix = cv2.imread('boat.png',0)
plt.imshow(im)
print("matrix shape: ",imMatrix.shape)
print("matrix chart:\n",imMatrix)
matrix shape: (512, 512) matrix chart: [[127 123 125 ... 165 169 166] [128 126 128 ... 169 163 167] [128 124 128 ... 178 160 175] ... [112 112 115 ... 101 97 104] [110 112 117 ... 104 93 105] [113 115 121 ... 102 95 97]]
In [11]:
#
imMatrixFloat = np.array(imMatrix,dtype="float")
U2, S2, V2 = svds(imMatrixFloat,k=150)
#U3, S3, V3 = svd(imMatrixFloat)
S2 = np.flip(S2)
print("U matrix:", U2)
print("V matrix:", V2)
print("singular values :", S2)
# print("U matrix:", U3)
# print("V matrix:", V3)
# print("singular values :", S3)
U matrix: [[-0.04314208 -0.0388528 -0.00897191 ... -0.00538612 0.07128615 -0.04873523] [-0.06252526 -0.02922972 -0.03186397 ... -0.00545031 0.07163346 -0.04872954] [-0.03363917 -0.01639127 -0.02515783 ... -0.0051772 0.07137113 -0.0487061 ] ... [-0.00683852 0.04003784 -0.00698873 ... -0.01915028 0.02201401 -0.04588184] [-0.05088047 0.08080697 0.01360937 ... -0.01518158 0.02571904 -0.04489299] [-0.0664585 0.02636251 0.00025653 ... -0.01133331 0.02647854 -0.04434909]] V matrix: [[-0.10177899 -0.02230333 -0.03907341 ... -0.03041493 0.01349579 0.06670612] [-0.01336242 -0.06118184 -0.0463375 ... 0.01118427 -0.01978009 -0.01719997] [ 0.02453641 0.0140836 -0.03390233 ... -0.02548265 0.071176 0.04830984] ... [ 0.00852677 0.00336625 -0.00123133 ... 0.02419142 0.02411709 0.02163343] [-0.03034815 -0.03016858 -0.02397761 ... 0.06212641 0.06084498 0.05969371] [-0.03914198 -0.04403926 -0.04302223 ... -0.03614463 -0.0344949 -0.03343196]] singular values : [68183.03140563 8317.75865645 7022.67647098 4965.23118021 4843.77635891 4672.11229171 3893.80139774 3763.67386444 2938.06976794 2687.95465976 2479.10387043 2395.69027861 2283.94361877 2057.77127137 1971.33864174 1754.25957317 1737.65852846 1694.5856797 1589.45682714 1567.50279264 1549.33418594 1462.83349689 1407.70395578 1369.63146203 1365.80840481 1308.46158837 1240.53162007 1236.93604373 1199.07757498 1142.86685071 1134.24183349 1110.27297279 1076.64071901 1064.65856912 1036.11482811 1000.25653759 987.17865842 968.39189308 943.65733303 938.36414711 919.81917896 912.66653824 896.47742752 882.22164959 866.10306276 856.32283548 837.92631 831.34777715 817.37554291 795.34371777 772.86856891 767.63689455 757.18993596 744.64223946 729.91839672 723.05569909 705.6042133 683.01465234 678.07861311 669.9984535 661.35082187 652.4959735 643.97856885 624.8021051 618.20710155 608.60781107 600.06333691 589.2780472 573.77669566 568.13624366 564.97963306 555.24559059 550.18333743 544.84562067 538.4338961 536.60309963 522.31762481 520.39756128 512.45204135 503.65130046 495.28682669 484.9982134 476.01704977 473.26592387 471.27243015 460.60482023 456.19717702 449.83009049 446.53473846 437.39401071 435.44592515 427.485253 425.63035915 416.03704273 407.84686126 407.39886843 404.11608999 396.78166612 394.77038543 388.58982097 384.81324061 381.47345821 379.36112684 375.78848033 370.63869285 369.31134702 366.26732583 360.02227111 356.89778626 352.76812803 346.54814415 344.51072748 343.31254011 333.23236205 331.65963366 330.02210245 325.25086458 322.41275805 320.41629136 315.10662014 309.88523555 307.98585249 305.04631466 302.67437858 301.14475468 297.0215458 294.04005509 289.43195433 286.84204491 285.69458981 282.51536578 281.16284223 279.32107438 275.26490595 269.68870857 268.64144238 265.50113266 260.66629761 258.06382111 256.87939207 256.07543014 251.94355266 249.27635586 246.46075024 244.97475587 243.62956371 241.60104126 239.58917507 235.10778471 234.14054425]
In [12]:
S150 = []
for i in range(150): #0:150 of singular values
S150.append(S2[i])
#print(S150)
In [13]:
#2a
S150_N = range(150)
#print(S150_N)
plt.plot(S150_N,S150)
plt.xlabel("N");
plt.ylabel("Singular Values");
plt.title("First 150 Singular Values from pSVD");
Regarding the number of terms needed to approximate the image, the singular values imply that it would be computationally cost-efficient and still accurate if we used less than 50 terms, as the singular values decrease as terms increase. The SVD function prints the singular values in decreasing order, while svds prints them in increasing order. Scipy documentation states that the svds function does not guarantee the order of the output of singular values upon computation of matrix. To correct this, I had to flip the values.
In [14]:
#b.
U8, S8, V8 = svds(imMatrixFloat,k=8)
#print(U8.shape,S8.shape,V8.shape)
m = len(U8[0])
n = len(S8)
S8Boat = np.zeros((m,n))
for i in range(n):
S8Boat[i,i] = S8[i]
S8Boat = U8@S8Boat@V8
plt.subplot(1,2,1);
plt.imshow(S8Boat, cmap="gray");
plt.title("Boat S=8");
plt.subplot(1,2,2);
plt.imshow(im);
plt.title("Original Boat");
#print("U matrix:", U8)
#print("V matrix:", V8)
#print("singular values :", S8)
In [15]:
U16, S16, V16 = svds(imMatrixFloat,k=16)
m = len(U16[0])
n = len(S16)
S16Boat = np.zeros((m,n))
for i in range(n):
S16Boat[i,i] = S16[i]
S16Boat = U16@S16Boat@V16
plt.subplot(1,2,1);
plt.imshow(S16Boat, cmap="gray");
plt.title("Boat S=16");
plt.subplot(1,2,2);
plt.imshow(im);
plt.title("Original Boat");
In [16]:
U32, S32, V32 = svds(imMatrixFloat,k=32)
m = len(U32[0])
n = len(S32)
S32Boat = np.zeros((m,n))
for i in range(n):
S32Boat[i,i] = S32[i]
S32Boat = U32@S32Boat@V32
plt.subplot(1,2,1);
plt.imshow(S32Boat, cmap="gray");
plt.title("Boat S=32");
plt.subplot(1,2,2);
plt.imshow(im);
plt.title("Original Boat");
In [17]:
U64, S64, V64 = svds(imMatrixFloat,k=64)
m = len(U64[0])
n = len(S64)
S64Boat = np.zeros((m,n))
for i in range(n):
S64Boat[i,i] = S64[i]
S64Boat = U64@S64Boat@V64
plt.subplot(1,2,1);
plt.imshow(S64Boat, cmap="gray");
plt.title("Boat S=64");
plt.subplot(1,2,2);
plt.imshow(im);
plt.title("Original Boat");
In [18]:
U128, S128, V128 = svds(imMatrixFloat,k=128)
m = len(U128[0])
n = len(S128)
S128Boat = np.zeros((m,n))
for i in range(n):
S128Boat[i,i] = S128[i]
S128Boat = U128@S128Boat@V128
plt.subplot(1,2,1);
plt.imshow(S128Boat, cmap="gray");
plt.title("Boat S=128");
plt.subplot(1,2,2);
plt.imshow(im);
plt.title("Original Boat");
The quality of the images increases substantially at each doubling of the rank. This makes sense, as the greater the rank, the higher the image quality should be. As I expected, the quality doesn't improve much from rank 50 on. Even at rank 8, without prior knowledge, it may be difficult for some to determine what the image is.
In [19]:
#2 C.
imMatrixBoat = plt.imread('boat.png')
MaxError = 0.005 #.5%
Error = 0
k = 1
A_k_Rel = 0
imMatrixFloat = np.array(imMatrixBoat,dtype="float")
while k < imMatrixFloat[0].size:
U, S, VT = svds(imMatrixFloat,k)
Ak = U@np.diag(S)@VT
A_k = np.linalg.norm(imMatrixFloat - Ak, ord="fro")
A_k2 = np.linalg.norm(imMatrixFloat, ord="fro")
A_k_Rel = A_k/A_k2
print(A_k_Rel,k)
k += 1
Error = A_k_Rel
if Error < MaxError:
break
0.25838308523418385 1 0.22994218704511368 2 0.20729973804088486 3 0.19499786097220406 4 0.1825220679177362 5 0.1700951654602372 6 0.16089986351273644 7 0.151806444619449 8 0.14598744544499206 9 0.14093242254289107 10 0.1364851332507304 11 0.13219705864980902 12 0.12817530234287106 13 0.12481537765934894 14 0.12165014106009751 15 0.11908394399055973 16 0.11651116503278984 17 0.11401050190004675 18 0.11176423159689004 19 0.1095354073555779 20 0.10731325752612714 21 0.10529277979874575 22 0.10338651521894042 23 0.10154900362664455 24 0.09968815408329011 25 0.09794918105342981 26 0.09635930306802148 27 0.09475218574005134 28 0.09321669042697209 29 0.09179952232373238 30 0.09038194419833898 31 0.08900246637118159 32 0.08768550054738738 33 0.08637826995759632 34 0.08512168558545839 35 0.08393363801844861 36 0.08276005932341608 37 0.08161479217910632 38 0.08051220437619587 39 0.07940689883907397 40 0.07833015828938342 41 0.07725543974894997 42 0.07620414479708364 43 0.07517200660395218 44 0.07416364680544955 45 0.07316450021514993 46 0.07219486469099255 47 0.0712275039216747 48 0.07027973054953261 49 0.06937042853321342 50 0.0685007118486341 51 0.06763177183393275 52 0.06677546998036984 53 0.0659367357858259 54 0.06512066665424107 55 0.06430980483205559 56 0.06352799173130723 57 0.0627866022904079 58 0.062047222295352075 59 0.061316757801219875 60 0.06059655876177062 61 0.059887198180105655 62 0.059188062985585294 63 0.058522314636191106 64 0.057863126754458344 65 0.057217004678829925 66 0.05658182547027866 67 0.05596244569492641 68 0.055368825692208475 69 0.054780574595652824 70 0.05419256282700216 71 0.05361851678308312 72 0.05304884698816271 73 0.05248417213211241 74 0.05192678285635042 75 0.05136719117270908 76 0.050831314688122985 77 0.05029372406788223 78 0.049766878632077224 79 0.04925262298542284 80 0.04875014726928096 81 0.04826341774671908 82 0.0477898602579015 83 0.04731710318390052 84 0.0468436094738252 85 0.04638679561537049 86 0.045934268853580436 87 0.045489969643349566 88 0.045047869906980946 89 0.044619567229933106 90 0.04419097546264614 91 0.04377394030916044 92 0.043356557185437475 93 0.04295398846183 94 0.04256352620716282 95 0.042170317967119064 96 0.04177980918124555 97 0.041399859666125687 98 0.041020286378605476 99 0.04064912455507789 100 0.04028182087654403 101 0.039917572641799806 102 0.039554048326035644 103 0.039194061447505396 104 0.03884067211054286 105 0.03848659936213011 106 0.03813513259231857 107 0.0377924444497478 108 0.03745262397200284 109 0.037117617634744776 110 0.036791428543352316 111 0.03646619708909653 112 0.03614032769826071 113 0.035830602088471754 114 0.035521130633773955 115 0.03521202787576743 116 0.03490917798572439 117 0.03460900896062075 118 0.03430996830771518 119 0.034018255958517374 120 0.03373373149131849 121 0.03345030835074078 122 0.033169916520726865 123 0.032891533275656255 124 0.032613616320536144 125 0.0323409660050787 126 0.03207151310612432 127 0.031808262456566316 128 0.03154756370926271 129 0.03128680027448944 130 0.031029689040131297 131 0.03077291628089567 132 0.030517378229598758 133 0.030267142938908257 134 0.030024981583957132 135 0.029782751236247445 136 0.029544233678121454 137 0.029312486962384136 138 0.029083552730347476 139 0.02885492371876773 140 0.028625914691994117 141 0.028402477610460265 142 0.02818203044836823 143 0.027964855496560085 144 0.027748622343730274 145 0.027533087054993695 146 0.027319467688351452 147 0.027107749517957126 148 0.026902302755168243 149 0.02669698139453377 150 0.026492152862080537 151 0.026290665752986952 152 0.02609323976308264 153 0.02589928433013123 154 0.025708889153416785 155 0.025518991121429857 156 0.02532862130197214 157 0.025142854623056987 158 0.02496056632561149 159 0.024778752417636742 160 0.024597561494659418 161 0.02441712913919912 162 0.024239178929089528 163 0.024066068634729022 164 0.02389227681705974 165 0.023720602476176928 166 0.0235491278535331 167 0.023382976937065984 168 0.02321677302442067 169 0.023051831406640834 170 0.022888095092896958 171 0.022724964915516923 172 0.02256358552425364 173 0.022404264675706803 174 0.022245089055115505 175 0.02208717024761301 176 0.021932460077517448 177 0.02178071727787504 178 0.021628565555693112 179 0.02148023556882659 180 0.02133313605874078 181 0.021185285031945025 182 0.02103814930962606 183 0.02089293364297572 184 0.020749299288369147 185 0.020606306556711345 186 0.020463578644739546 187 0.020321005032139447 188 0.02017986414577341 189 0.020039779194144416 190 0.01990193635455364 191 0.01976491100087897 192 0.019629003079415893 193 0.01949395242373483 194 0.01935921607555706 195 0.019225933344575893 196 0.019093645433502317 197 0.018961480834593578 198 0.01882890248818856 199 0.018696628372210488 200 0.018565751261378622 201 0.018436423092543005 202 0.018308798728858443 203 0.018182397873909725 204 0.018056891689406968 205 0.017931118994970595 206 0.017805721554111125 207 0.01768048735500052 208 0.017556684548460984 209 0.017434499398222246 210 0.017313078715290992 211 0.01719253429501144 212 0.01707267015443439 213 0.01695305490714494 214 0.016837446362020277 215 0.0167219293979673 216 0.016605744055655125 217 0.01649021665564817 218 0.016375493399773136 219 0.01626044318786747 220 0.01614604620780916 221 0.01603247690630429 222 0.01591902705339814 223 0.01580741062742309 224 0.015696549623075456 225 0.015585871091767125 226 0.015474693057856466 227 0.015364652406941138 228 0.015255441352938565 229 0.015148126714388532 230 0.015042278046284363 231 0.014936320831368437 232 0.01483098258488186 233 0.014725846891013071 234 0.014620522630251917 235 0.014516711761363089 236 0.014412854153370156 237 0.014311213252237201 238 0.014209747411686654 239 0.014109050992967786 240 0.014008501426603011 241 0.013908669646135624 242 0.013809534716751488 243 0.013711690409443223 244 0.013614465101782448 245 0.013517144727737556 246 0.013421132180228182 247 0.013324987717915019 248 0.0132291844500209 249 0.01313330920502808 250 0.013037836050944369 251 0.01294425288313353 252 0.012850477358268947 253 0.012756887092844173 254 0.012663459845527535 255 0.012570905060496124 256 0.012478389180098817 257 0.012387158285214953 258 0.012297097111779308 259 0.01220778168808225 260 0.012118688061964856 261 0.012029834915367138 262 0.011940467600451242 263 0.0118514996532998 264 0.011763123827237612 265 0.011674970513477415 266 0.011588175019932681 267 0.011501350840288572 268 0.0114148869924179 269 0.011329404614903775 270 0.011244075435309345 271 0.011158924795712684 272 0.01107422748943704 273 0.010989087770940547 274 0.010904666209291132 275 0.010820101486413028 276 0.01073609634320853 277 0.010652787712224412 278 0.010570336861079718 279 0.010488514180315442 280 0.010406798394534294 281 0.010325629807259408 282 0.01024427305518904 283 0.010162765585560361 284 0.010083604801184719 285 0.010004231991563138 286 0.009925018923103634 287 0.009846175158744811 288 0.009767261327807673 289 0.009688625093112477 290 0.009609695274107257 291 0.009531452422080835 292 0.009452876630441892 293 0.00937580909522363 294 0.009298991369346733 295 0.009222769948593596 296 0.009147319370967728 297 0.009072195817811456 298 0.008997824221156938 299 0.008923205010703327 300 0.008849175691033882 301 0.00877491822656333 302 0.008700506738005369 303 0.008627118860460047 304 0.008553914874583785 305 0.008481220426069943 306 0.00840829747094272 307 0.008337032767279408 308 0.008265977914642961 309 0.008195199846549057 310 0.008124701371581652 311 0.008054162736646978 312 0.00798465095036982 313 0.007916180142765784 314 0.007847170612257752 315 0.007778139755700048 316 0.007709186966292789 317 0.007641104989594354 318 0.007573063110835177 319 0.007505068930429662 320 0.007437209164985939 321 0.00736935603085766 322 0.0073024471146105295 323 0.007236357843680375 324 0.007170023076686261 325 0.007104074867927581 326 0.007038472234531549 327 0.006973687666000785 328 0.006908580059377432 329 0.006843642903053732 330 0.006779117601029213 331 0.006714424525843982 332 0.006650435092981372 333 0.006587202736727909 334 0.006524844077348323 335 0.006462015126814522 336 0.006399577488628601 337 0.006337343446187977 338 0.006275561350410215 339 0.0062139550217604495 340 0.00615229223703255 341 0.00609084320628955 342 0.006029831288115075 343 0.005968946394718497 344 0.005908635525308089 345 0.005848156579810147 346 0.005788140243025763 347 0.005729849695292073 348 0.00567141405346327 349 0.0056129491102703466 350 0.005556283130906936 351 0.005499470998787884 352 0.005442329780182792 353 0.005385431877328408 354 0.005329137543707406 355 0.0052731668335453665 356 0.005217412253178543 357 0.005161290169037297 358 0.005105159857752261 359 0.005049849596372872 360 0.00499475851622371 361
2c. The storage required for this data compression is actually more than it would take for the image itself due to the large rank. For scenarios when the rank is smaller, the storage would likely be more efficient. According to low-rank approximations notes, we can rewrite the SVD of A as s_i * u_i * (v_i).Transpose, where u_i is the ith column of u, v_i is the ith column of v, and s_i is the ith number in s. Based on this, we can calculate the computational costs.
In [20]:
print("svd cost: ",361*512 + 361*512 + 361)
print("non-svd image cost: ",512*512)
svd cost: 370025 non-svd image cost: 262144
In [44]:
3.
Mandrill = plt.imread('mandrill.png',3)
#plt.imshow(Mandrill)
#imMatrixFloat = np.array(Mandrill,dtype="float")
#print("matrix chart:\n",imMatrixFloat)
B = Mandrill[:,:,0]
G = Mandrill[:,:,1]
R = Mandrill[:,:,2]
MandrillMerged = cv2.merge((B,G,R))
#print("B\n",B);print("G\n",G);print("R\n",R)
#print(imMatrixFloat.shape)
#print(B.shape)
B = np.array(B,dtype="float")
G = np.array(G,dtype="float")
R = np.array(R,dtype="float")
plt.subplot(1,4,1);
plt.imshow(MandrillMerged);
plt.axis("off");
plt.subplot(1,4,2);
plt.imshow(B,cmap="Blues");
plt.axis("off");
plt.subplot(1,4,3);
plt.imshow(G,cmap="Greens");
plt.axis("off");
plt.subplot(1,4,4);
plt.imshow(R,cmap="Reds");
plt.axis("off");
In [46]:
U, SB, V = svds(B,k=150)
U, SG, V = svds(G,k=150)
U, SR, V = svds(R,k=150)
SB = np.flip(SB)
SG = np.flip(SG)
SR = np.flip(SR)
print("Blue singular values :", SB)
print("Green singular values :", SG)
print("Red singular values :", SR)
Blue singular values : [286.13805986 43.77857877 25.94407897 23.21779687 20.42378224 15.36809613 11.70547229 9.50867542 8.67742002 8.52962908 8.06470169 7.92829338 7.49241713 7.21256017 7.07859225 6.80481974 6.78702554 6.59775002 6.34444912 6.27000469 6.08477848 6.0651322 5.98112968 5.90912282 5.79641223 5.71631769 5.55400673 5.48549741 5.4278859 5.36949971 5.23865518 5.21449686 5.15996747 5.0978847 5.06343942 4.98774865 4.96490783 4.91728274 4.87225009 4.78972004 4.69754733 4.66516373 4.63223619 4.59160083 4.53199128 4.48984522 4.44283545 4.4047143 4.32421821 4.30008706 4.27935147 4.23873296 4.18001098 4.13828227 4.09541386 4.07143158 4.0161852 3.9948373 3.97596539 3.92044294 3.89176996 3.84634479 3.8339429 3.76078728 3.73634556 3.69767351 3.6906152 3.62851114 3.6097143 3.57746301 3.54680299 3.5288602 3.45921427 3.4536668 3.4398161 3.4192336 3.35930297 3.35148371 3.34055154 3.3131483 3.2703594 3.24927525 3.24159485 3.22054546 3.16737474 3.16384335 3.14141357 3.11973828 3.07743748 3.05507506 3.03104943 3.02909626 3.00205105 2.98400013 2.9682056 2.94581549 2.91802381 2.90625514 2.85753358 2.82491373 2.80626318 2.7974372 2.7795233 2.76641676 2.73494614 2.71700835 2.68535128 2.67892627 2.65653035 2.61637559 2.61372918 2.58531146 2.58363705 2.56123323 2.54532764 2.52326978 2.49173892 2.4692303 2.46455078 2.4316565 2.41522727 2.4027983 2.39554944 2.38705176 2.37222111 2.34606743 2.32803234 2.30532967 2.30477068 2.28113027 2.26957189 2.25112235 2.2414308 2.22544111 2.2004578 2.19965868 2.17814127 2.17372355 2.142466 2.13072641 2.12125392 2.08429473 2.07687699 2.06542805 2.04603846 2.03191058 2.02107823 2.0194466 2.01157926 2.0018272 ] Green singular values : [262.49823941 43.85032798 24.71768069 19.87212148 16.63209544 14.37883925 12.13199915 10.47151231 10.37152122 9.49684037 9.10673572 8.64389425 7.98869037 7.95108445 7.89844725 7.742241 7.51300828 7.32293911 7.22259141 7.11540347 6.90825934 6.89464841 6.71322576 6.66745633 6.58998969 6.52014775 6.47967722 6.43215266 6.34637331 6.27182828 6.23626199 6.15722682 6.04785083 5.99964527 5.88309991 5.79725618 5.7237147 5.70086162 5.68758985 5.58091491 5.5169807 5.50775469 5.3783189 5.3518015 5.32522965 5.2766939 5.22722341 5.16302273 5.11829114 5.08607244 5.05679392 5.02886479 4.94379065 4.92550721 4.88779808 4.80203362 4.77257634 4.75571326 4.6909296 4.64728199 4.61100152 4.56890732 4.53003106 4.49591174 4.45686483 4.41304032 4.38309281 4.33486422 4.3196131 4.27482673 4.25605515 4.2391709 4.22631427 4.16009742 4.13999098 4.11612348 4.09297045 4.02225711 4.00928424 3.95362 3.93991379 3.92766021 3.915017 3.84864559 3.8126066 3.80786351 3.77628158 3.74979701 3.73867964 3.70353885 3.64851459 3.64006304 3.59877617 3.59000454 3.56923341 3.53272269 3.49640019 3.46395304 3.43893753 3.43217909 3.40080912 3.37686546 3.34962235 3.31830739 3.30625683 3.30119327 3.25930489 3.22827826 3.2257985 3.1815791 3.17196462 3.13853565 3.11398527 3.10764931 3.04358946 3.0249328 3.00775234 3.00033529 2.97278926 2.95083986 2.91352177 2.89954022 2.88429269 2.86301168 2.85005267 2.82667008 2.81163458 2.7940243 2.77041301 2.75375777 2.74289252 2.7296874 2.7074844 2.69689429 2.67953108 2.66271945 2.63195914 2.59994788 2.58441879 2.57592847 2.5561297 2.55129363 2.53177354 2.51498602 2.47503003 2.45028606 2.43778112 2.42591622 2.40200917 2.38078357] Red singular values : [240.64085155 52.12753585 30.75756124 23.04534025 18.23008113 17.07265284 15.34383014 13.12058849 11.40070621 11.23306226 10.89875376 10.31678195 9.63014546 8.9311381 8.81654578 8.41753194 8.23314809 7.83137162 7.54221018 7.51869863 7.20573899 7.11117285 6.99633605 6.91100615 6.80264495 6.68924383 6.59690757 6.51694186 6.40361371 6.2765096 6.17717885 6.13453545 6.06895528 5.91972766 5.85354446 5.79110712 5.72620891 5.65918773 5.60415144 5.57723415 5.53114603 5.4614575 5.39386218 5.36301436 5.30182457 5.26945006 5.23558588 5.17520239 5.08419561 5.0755036 5.06582308 4.97834715 4.91622693 4.90195264 4.82442364 4.73435004 4.72566592 4.66747627 4.6631467 4.62371194 4.6140662 4.54627762 4.50006531 4.49374405 4.46176875 4.42500792 4.41872482 4.31784992 4.25819875 4.23557562 4.23221877 4.18117237 4.16886205 4.11781476 4.08112887 4.04507294 4.01455836 3.95934708 3.94252671 3.90870909 3.88531789 3.84357057 3.80759212 3.77075415 3.73983186 3.72908495 3.69413534 3.66060164 3.62442138 3.60064893 3.58006202 3.57782314 3.53029582 3.50239322 3.47306503 3.45893761 3.4313387 3.4098006 3.3986317 3.38914381 3.34450461 3.31162694 3.30761503 3.25570496 3.23920737 3.22720382 3.20612585 3.17363392 3.13572212 3.11282101 3.0897685 3.06614276 3.03393697 3.02396977 3.0080281 3.00092857 2.98451005 2.96363247 2.91757251 2.8997601 2.87073033 2.85860888 2.83075172 2.82789054 2.82396343 2.79305496 2.78123888 2.74998483 2.73548817 2.72996891 2.71784336 2.69337384 2.6780351 2.65998441 2.64239633 2.62511529 2.60820706 2.60214036 2.56267537 2.55695183 2.53500169 2.52519976 2.5073102 2.48376491 2.46795971 2.45564485 2.44511596 2.42486778 2.41953883 2.38260879]
The singular values for red, green, and blue decay roughly as fast as each other. As seen above, blue has the largest initial singular value, but also decays the fastest, as the 150th singular value for blue is the lowest singular value.
In [62]:
# Convert to uint8 in [0,255] if needed
if Mandrill.max() <= 1.0:
Mandrill_uint8 = (Mandrill * 255).round().astype(np.uint8)
else:
Mandrill_uint8 = Mandrill.astype(np.uint8)
R = Mandrill_uint8[:, :, 0]
G = Mandrill_uint8[:, :, 1]
B = Mandrill_uint8[:, :, 2]
def low_rank_svd_dense(channel, k):
U, S, Vt = svd(channel, full_matrices=False)
U_k = U[:, :k]
S_k = S[:k]
Vt_k = Vt[:k, :]
S_matrix = np.diag(S_k)
approx = U_k @ S_matrix @ Vt_k
approx = np.clip(approx, 0, 255)
return approx.astype(np.uint8)
k = 8
Red8 = low_rank_svd_dense(R, k)
Green8 = low_rank_svd_dense(G, k)
Blue8 = low_rank_svd_dense(B, k)
Mandrill_approx = np.stack((Red8, Green8, Blue8), axis=2)
plt.figure(figsize=(6, 3))
plt.subplot(1, 2, 1)
plt.imshow(Mandrill_approx)
plt.title(f"Mandrill S={k}")
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(Mandrill_uint8)
plt.title("Original Mandrill")
plt.axis('off')
plt.tight_layout()
plt.show()
In [63]:
# Ensure uint8, [0,255]
if Mandrill.max() <= 1.0:
Mandrill_uint8 = (Mandrill * 255).round().astype(np.uint8)
else:
Mandrill_uint8 = Mandrill.astype(np.uint8)
R = Mandrill_uint8[:, :, 0]
G = Mandrill_uint8[:, :, 1]
B = Mandrill_uint8[:, :, 2]
def low_rank_svd_dense(channel, k):
U, S, Vt = svd(channel, full_matrices=False)
U_k = U[:, :k]
S_k = S[:k]
Vt_k = Vt[:k, :]
S_matrix = np.diag(S_k)
approx = U_k @ S_matrix @ Vt_k
approx = np.clip(approx, 0, 255)
return approx.astype(np.uint8)
k = 16
Red16 = low_rank_svd_dense(R, k)
Green16 = low_rank_svd_dense(G, k)
Blue16 = low_rank_svd_dense(B, k)
Mandrill_approx = np.stack((Red16, Green16, Blue16), axis=2)
plt.figure(figsize=(6, 3))
plt.subplot(1, 2, 1)
plt.imshow(Mandrill_approx)
plt.title(f"Mandrill S={k}")
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(Mandrill_uint8)
plt.title("Original Mandrill")
plt.axis('off')
plt.tight_layout()
plt.show()
In [65]:
R = Mandrill_uint8[:, :, 0]
G = Mandrill_uint8[:, :, 1]
B = Mandrill_uint8[:, :, 2]
def low_rank_svd_dense(channel, k):
U, S, Vt = svd(channel, full_matrices=False)
U_k = U[:, :k]
S_k = S[:k]
Vt_k = Vt[:k, :]
S_matrix = np.diag(S_k)
approx = U_k @ S_matrix @ Vt_k
approx = np.clip(approx, 0, 255)
return approx.astype(np.uint8)
k = 32
Red32 = low_rank_svd_dense(R, k)
Green32 = low_rank_svd_dense(G, k)
Blue32 = low_rank_svd_dense(B, k)
Mandrill_approx = np.stack((Red32, Green32, Blue32), axis=2)
plt.figure(figsize=(6, 3))
plt.subplot(1, 2, 1)
plt.imshow(Mandrill_approx)
plt.title(f"Mandrill S={k}")
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(Mandrill_uint8)
plt.title("Original Mandrill")
plt.axis('off')
plt.tight_layout()
plt.show()
In [66]:
R = Mandrill_uint8[:, :, 0]
G = Mandrill_uint8[:, :, 1]
B = Mandrill_uint8[:, :, 2]
def low_rank_svd_dense(channel, k):
U, S, Vt = svd(channel, full_matrices=False)
U_k = U[:, :k]
S_k = S[:k]
Vt_k = Vt[:k, :]
S_matrix = np.diag(S_k)
approx = U_k @ S_matrix @ Vt_k
approx = np.clip(approx, 0, 255)
return approx.astype(np.uint8)
k = 64
Red64 = low_rank_svd_dense(R, k)
Green64 = low_rank_svd_dense(G, k)
Blue64 = low_rank_svd_dense(B, k)
Mandrill_approx = np.stack((Red64, Green64, Blue64), axis=2)
plt.figure(figsize=(6, 3))
plt.subplot(1, 2, 1)
plt.imshow(Mandrill_approx)
plt.title(f"Mandrill S={k}")
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(Mandrill_uint8)
plt.title("Original Mandrill")
plt.axis('off')
plt.tight_layout()
plt.show()
In [67]:
R = Mandrill_uint8[:, :, 0]
G = Mandrill_uint8[:, :, 1]
B = Mandrill_uint8[:, :, 2]
def low_rank_svd_dense(channel, k):
U, S, Vt = svd(channel, full_matrices=False)
U_k = U[:, :k]
S_k = S[:k]
Vt_k = Vt[:k, :]
S_matrix = np.diag(S_k)
approx = U_k @ S_matrix @ Vt_k
approx = np.clip(approx, 0, 255)
return approx.astype(np.uint8)
k = 128
Red128 = low_rank_svd_dense(R, k)
Green128 = low_rank_svd_dense(G, k)
Blue128 = low_rank_svd_dense(B, k)
Mandrill_approx = np.stack((Red128, Green128, Blue128), axis=2)
plt.figure(figsize=(6, 3))
plt.subplot(1, 2, 1)
plt.imshow(Mandrill_approx)
plt.title(f"Mandrill S={k}")
plt.axis('off')
plt.subplot(1, 2, 2)
plt.imshow(Mandrill_uint8)
plt.title("Original Mandrill")
plt.axis('off')
plt.tight_layout()
plt.show()
As the rank increases, the image becomes more accurate. This is expected.
In [ ]: