Jeremy Knight Math 639 2/22/11 Householder's Method To obtain a symmetric tridiagonal matrix LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkobWZlbmNlZEdGJDYkLUYjNictRi82JVEibkYnRjJGNS1JI21vR0YkNi1RIitGJy9GNlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkovJSlzdHJldGNoeUdGSi8lKnN5bW1ldHJpY0dGSi8lKGxhcmdlb3BHRkovJS5tb3ZhYmxlbGltaXRzR0ZKLyUnYWNjZW50R0ZKLyUnbHNwYWNlR1EsMC4yMjIyMjIyZW1GJy8lJ3JzcGFjZUdGWS1JI21uR0YkNiRRIjFGJ0ZGRjJGNUZGRjJGNS8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGRg== similar to the symmetric matrix LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JJW1zdXBHRiQ2JUYrLUYjNiUtSShtZmVuY2VkR0YkNiQtRiM2JS1JI21uR0YkNiRRIjFGJ0Y5Ri9GMkY5Ri9GMi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRictRjY2LVEifkYnRjlGO0Y+RkBGQkZERkZGSC9GS1EmMC4wZW1GJy9GTkZeb0Y5 construct the following matrices LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYtLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkobWZlbmNlZEdGJDYkLUYjNiUtSSNtbkdGJDYkUSIyRicvRjZRJ25vcm1hbEYnRjJGNUZDRjJGNS8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRictSSNtb0dGJDYtUSIsRidGQy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R0ZOLyUqc3ltbWV0cmljR0ZOLyUobGFyZ2VvcEdGTi8lLm1vdmFibGVsaW1pdHNHRk4vJSdhY2NlbnRHRk4vJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR1EsMC4zMzMzMzMzZW1GJy1GLDYlRi4tRiM2JS1GOzYkLUYjNiUtRkA2JFEiM0YnRkNGMkY1RkNGMkY1RkVGSC1GLDYlRi4tRiM2JS1GOzYkLUYjNiUtRkA2JFEiNEYnRkNGMkY1RkNGMkY1RkVGSC1GSTYtUSMuLkYnRkNGTC9GUEZORlFGU0ZVRldGWS9GZm5RLDAuMjIyMjIyMmVtRicvRmluRmduLUZJNi1RIi5GJ0ZDRkxGZHBGUUZTRlVGV0ZZRmVuRmdwRkgtRiw2JUYuLUYjNiUtRjs2JC1GIzYnLUYvNiVRIm5GJ0YyRjUtRkk2LVEoJm1pbnVzO0YnRkNGTEZkcEZRRlNGVUZXRllGZXAvRmluRmZwLUZANiRRIjFGJ0ZDRjJGNUZDRjJGNUZFRkM= , where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1cEdGJDYlLUkjbWlHRiQ2JVEiQUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkobWZlbmNlZEdGJDYkLUYjNiUtRi82JVEia0YnRjJGNUYyRjUvRjZRJ25vcm1hbEYnRjJGNS8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRictSSNtb0dGJDYtUSI9RidGQi8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGTS8lKXN0cmV0Y2h5R0ZNLyUqc3ltbWV0cmljR0ZNLyUobGFyZ2VvcEdGTS8lLm1vdmFibGVsaW1pdHNHRk0vJSdhY2NlbnRHRk0vJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0Zmbi1GOzYkLUYjNiQtSShtc3Vic3VwR0YkNigtRi82JVEiYUYnRjJGNS1GIzYlLUYvNiVRI2lqRidGMkY1RjJGNUY4RkQvJS9zdWJzY3JpcHRzaGlmdEdGRi9JK21zZW1hbnRpY3NHRiRRLFtub25lLG5vbmVdRidGQkZCRkI= for each LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYvLUkjbWlHRiQ2JVEia0YnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JI21uR0YkNiRRIjFGJ0Y5LUY2Ni1RIixGJ0Y5RjsvRj9GMUZARkJGREZGRkgvRktRJjAuMGVtRicvRk5RLDAuMzMzMzMzM2VtRictRlA2JFEiMkYnRjlGUy1GNjYtUSMuLkYnRjlGO0Y+RkBGQkZERkZGSC9GS1EsMC4yMjIyMjIyZW1GJy9GTkZYLUY2Ni1RIi5GJ0Y5RjtGPkZARkJGREZGRkhGV0Zdb0ZTLUYsNiVRIm5GJ0YvRjItRjY2LVEoJm1pbnVzO0YnRjlGO0Y+RkBGQkZERkZGSEZbby9GTkZcb0ZPRjk=.
<Text-field style="Heading 1" layout="Heading 1">Householder(n,A0,display)</Text-field> Householder:=proc(n,A0,display) #n=dimension, A=nxn symmetric matrix local A,i,j,k,l,q,u,v,w,z,alpha,PROD,temp,RSQ: #option trace: A:=A0: for k from 1 to n-2 do q:=0: for j from (k+1) to n do q:=q+A[j,k]^2: end do: if A[k+1,k]=0 then alpha:=-q^(1/2); else alpha:=-(q^(1/2)*A[k+1,k])/(abs(A[k+1,k])); end if; #RSQ=2r^2 RSQ:=alpha^2-alpha*A[k+1,k]: #w=(1/sqrt(1RSQ))*v = 1/2r * v v:=Array(1..n): v[k+1]:=A[k+1,k]-alpha: for j from k+2 to n do v[j]:=A[j,k]: end do: #u=(1/RSQ)*A*v = (1/2r^2)*A*v = 1/r*A*w u:=Array(1..n): for j from k to n do temp:=0: for i from k+1 to n do temp:=temp+A[j,i]*v[i]: end do: u[j]:=(1/RSQ)*temp: end do: #PROD = v*u = (1/2r^2)*v*A*v = (1/r)Aw PROD:=0: for i from k+1 to n do PROD:=PROD+v[i]*u[i]: end do: #z = u-(1/2RSQ)vuv = u-(1/4r^2)vuv = u-wwu = (1/r)Aw-ww(1/r)Aw z:=Array(1..n): for j from k to n do z[j]:=u[j]-(PROD/(2*RSQ))*v[j]: end do; #A^(k+1) = A^(k)-vz-zv=(I-2ww)A^(k)(I-2ww) #Note: at each step the values of A are overwritten to create A^(k+1) for l from k+1 to n-1 do for j from l+1 to n do A[j,l]:=A[j,l]-v[l]*z[j]-v[j]*z[l]: A[l,j]:=A[j,l]: end do: A[l,l]:=A[l,l]-2*v[l]*z[l]: end do: A[n,n]:=A[n,n]-2*v[n]*z[n]: for j from k+2 to n do A[k,j]:=0: A[j,k]:=0: end do: A[k+1,k]:=A[k+1,k]-v[k+1]*z[k]: A[k,k+1]:=A[k+1,k]: #All other elements of A^(k+1) are the same as A^(k) if display=1 then print(k,A): end if: end do: return(A): end proc:
LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=
<Text-field style="Heading 1" layout="Heading 1">Householder2(n,A0,display)</Text-field> LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYvLUkjbWlHRiQ2JVElVGhpc0YnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYtUSJ+RidGMi8lJmZlbmNlR0YxLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRjEvJSpzeW1tZXRyaWNHRjEvJShsYXJnZW9wR0YxLyUubW92YWJsZWxpbWl0c0dGMS8lJ2FjY2VudEdGMS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRkktRiw2JVEodmVyc2lvbkYnRi9GMkY1LUYsNiVRJXdpbGxGJ0YvRjJGNS1GLDYlUSZwcmludEYnRi9GMkY1LUYsNiVRJHRoZUYnRi9GMkY1LUYsNiVRKW1hdHJpY2VzRidGL0YyRjVGMg==LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYrLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiUEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RIixGJ0Y+LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y0LyUpc3RyZXRjaHlHRkkvJSpzeW1tZXRyaWNHRkkvJShsYXJnZW9wR0ZJLyUubW92YWJsZWxpbWl0c0dGSS8lJ2FjY2VudEdGSS8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHUSwwLjMzMzMzMzNlbUYnLUYsNiVGLi1GIzYlLUY7NiRRIjJGJ0Y+RjJGNUZARkMtRiw2JUYuLUYjNiUtRjs2JFEiM0YnRj5GMkY1RkBGQy1GRDYtUSMuLkYnRj5GRy9GS0ZJRkxGTkZQRlJGVC9GV1EsMC4yMjIyMjIyZW1GJy9GWkZYLUZENi1RIi5GJ0Y+RkdGZ29GTEZORlBGUkZURlZGam9GPg== where LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYtLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEiUEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0YyRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RJyZzZG90O0YnRj4vJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkkvJSlzdHJldGNoeUdGSS8lKnN5bW1ldHJpY0dGSS8lKGxhcmdlb3BHRkkvJS5tb3ZhYmxlbGltaXRzR0ZJLyUnYWNjZW50R0ZJLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGWC1GLDYlRi4tRiM2JS1GOzYkUSIyRidGPkYyRjVGQEZDLUZENi1RIy4uRidGPkZHRkpGTEZORlBGUkZUL0ZXUSwwLjIyMjIyMjJlbUYnRlktRkQ2LVEiLkYnRj5GR0ZKRkxGTkZQRlJGVEZWRllGQy1GLDYlRi4tRiM2JS1GLzYlUSJuRidGMkY1RjJGNUZALUZENi1RIj1GJ0Y+RkdGSkZMRk5GUEZSRlQvRldRLDAuMjc3Nzc3OGVtRicvRlpGX3AtRiw2JS1GLzYlUSJRRidGMkY1RmZvRkBGPg== Householder2:=proc(n,A0,display) #n=dimension, A=nxn symmetric matrix #option trace; local A,i,j,k,l,P,q,u,v,w,z,alpha,PROD,temp,RSQ; A:=A0: for k from 1 to n-2 do q:=0: for j from (k+1) to n do q:=q+A[j,k]^2: end do: if A[k+1,k]=0 then alpha:=-q^(1/2); else alpha:=-(q^(1/2)*A[k+1,k])/(abs(A[k+1,k])); end if; #RSQ=2r^2 RSQ:=alpha^2-alpha*A[k+1,k]: #w=(1/sqrt(RSQ))*v = 1/2r * v v:=Array(1..n): w:=Array(1..n): v[k+1]:=A[k+1,k]-alpha: w[k+1]:=v[k+1]/(2*(0.5*RSQ)^(1/2)); for j from k+2 to n do v[j]:=A[j,k]: w[j]:=v[j]/(2*(0.5*RSQ)^(1/2)): end do: #Compute P^(k)=I-2w*w P:=Array(1..n,1..n): for i from 1 to n do for j from 1 to n do if i=j then P[i,j]:=1-2*w[i]*w[j]: else P[i,j]:=0-2*w[i]*w[j]: end if: end do: end do: print("v",v); print("w",w); print("P^",k,P); #u=(1/RSQ)*A*v = (1/2r^2)*A*v = 1/r*A*w u:=Array(1..n): for j from k to n do temp:=0: for i from k+1 to n do temp:=temp+A[j,i]*v[i]: end do: u[j]:=(1/RSQ)*temp: end do: #PROD = v*u = (1/2r^2)*v*A*v = (1/r)Aw PROD:=0: for i from k+1 to n do PROD:=PROD+v[i]*u[i]: end do: #z = u-(1/2RSQ)vuv = u-(1/4r^2)vuv = u-wwu = (1/r)Aw-ww(1/r)Aw z:=Array(1..n): for j from k to n do z[j]:=u[j]-(PROD/(2*RSQ))*v[j]: end do; #A^(k+1) = A^(k)-vz-zv=(I-2ww)A^(k)(I-2ww) #Note: at each step the values of A are overwritten to create A^(k+1) for l from k+1 to n-1 do for j from l+1 to n do A[j,l]:=A[j,l]-v[l]*z[j]-v[j]*z[l]: A[l,j]:=A[j,l]: end do: A[l,l]:=A[l,l]-2*v[l]*z[l]: end do: A[n,n]:=A[n,n]-2*v[n]*z[n]: for j from k+2 to n do A[k,j]:=0: A[j,k]:=0: end do: A[k+1,k]:=A[k+1,k]-v[k+1]*z[k]: A[k,k+1]:=A[k+1,k]: #All other elements of A^(k+1) are the same as A^(k) if display=1 then print(k,A): end if: end do: return(A); end proc:
JSFH
<Text-field style="Heading 1" layout="Heading 1">Householder3(n,A0,display)</Text-field> This version will print the values of T and Q such that LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2JVEiVEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIj1GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGTC1JJW1zdXBHRiQ2JS1GLDYlUSRRQVFGJ0YvRjItRiM2JS1GLDYlUSJ0RidGL0YyRi9GMi8lMXN1cGVyc2NyaXB0c2hpZnRHUSIwRidGOQ== Householder3:=proc(n,A0,display) #n=dimension, A=nxn symmetric matrix #option trace; local A,i,j,k,l,P,q,Q,u,v,w,z,alpha,PROD,temp,RSQ,SqrMatrixProduct; SqrMatrixProduct:=proc(n,A,B) #multiplies two square nxn matrices A and B local i,j,k,C: C:=Array(1..n,1..n): for i from 1 to n do for j from 1 to n do for k from 1 to n do C[i,k]:=C[i,k]+A[i,j]*B[j,k]; end do: end do: end do: return(C); end proc: Q:=Array(1..n,1..n): A:=A0: for k from 1 to n-2 do q:=0: for j from (k+1) to n do q:=q+A[j,k]^2: end do: if A[k+1,k]=0 then alpha:=-q^(1/2); else alpha:=-(q^(1/2)*A[k+1,k])/(abs(A[k+1,k])); end if; #RSQ=2r^2 RSQ:=alpha^2-alpha*A[k+1,k]: #w=(1/sqrt(RSQ))*v = 1/2r * v v:=Array(1..n): w:=Array(1..n): v[k+1]:=A[k+1,k]-alpha: w[k+1]:=v[k+1]/(2*(0.5*RSQ)^(1/2)); for j from k+2 to n do v[j]:=A[j,k]: w[j]:=v[j]/(2*(0.5*RSQ)^(1/2)): end do: #Compute P^(k)=I-2w*w P:=Array(1..n,1..n): for i from 1 to n do for j from 1 to n do if i=j then P[i,j]:=1-2*w[i]*w[j]: else P[i,j]:=0-2*w[i]*w[j]: end if: end do:k end do: # print("v",v); # print("w",w); # print("P^",k,P); if k=1 then Q:=P: else Q:=SqrMatrixProduct(n,Q,P): end if: #u=(1/RSQ)*A*v = (1/2r^2)*A*v = 1/r*A*w u:=Array(1..n): for j from k to n do temp:=0: for i from k+1 to n do temp:=temp+A[j,i]*v[i]: end do: u[j]:=(1/RSQ)*temp: end do: #PROD = v*u = (1/2r^2)*v*A*v = (1/r)Aw PROD:=0: for i from k+1 to n do PROD:=PROD+v[i]*u[i]: end do: #z = u-(1/2RSQ)vuv = u-(1/4r^2)vuv = u-wwu = (1/r)Aw-ww(1/r)Aw z:=Array(1..n): for j from k to n do z[j]:=u[j]-(PROD/(2*RSQ))*v[j]: end do; #A^(k+1) = A^(k)-vz-zv=(I-2ww)A^(k)(I-2ww) #Note: at each step the values of A are overwritten to create A^(k+1) for l from k+1 to n-1 do for j from l+1 to n do A[j,l]:=A[j,l]-v[l]*z[j]-v[j]*z[l]: A[l,j]:=A[j,l]: end do: A[l,l]:=A[l,l]-2*v[l]*z[l]: end do: A[n,n]:=A[n,n]-2*v[n]*z[n]: for j from k+2 to n do A[k,j]:=0: A[j,k]:=0: end do: A[k+1,k]:=A[k+1,k]-v[k+1]*z[k]: A[k,k+1]:=A[k+1,k]: #All other elements of A^(k+1) are the same as A^(k) if display=1 then print(k,A): end if: end do: return(A,Q); end proc: