2.
What is ||A|| ?The following solution was supplied by Brian Medley
Siam Challenge Problem 3
(1) From the start, I assumed that this can be tackled by finding the
limit (as n->infinity) of the norms of the n*n matrices formed by
truncating A after the nth row and nth column, (matrices mude up of the top
left corner of A). If you don't wish to read my (non-rigorous) reasons for
trusting this assumption, skip straight to (2) below.
Using such a matrix (top corner of A) as a kind of "approximation" to A
would have an effect, when applied to an element of l_2, equivalent to the
following:
Truncate the element of l_2 after its nth postition, apply the matrix A to
this, and then truncate the result after its nth position.
I convinced myself that the two truncations would not matter two much if n
is large: the values removed from the later positions of original element
would be multiplied by the small numbers in A's columns n+1 and beyond, and
so make only small contributions to the result, whilst the removed values
of the image element would be small numbers anyway, having been generated
by calculations involving the small values in A's lower rows.
Thus an l_2- vector which is stretched a lot by A will be stretched a lot
if truncated and treated with a top-corner approximation to A, so the norms
of the top-corner sub-matrices should tend to the norm of A. Finally, the
l_2 norm is the infinite equivalent of the Euclidean norm on R^n, so in my
work on finite mtrices I can work with lengths of vectors in the usual way.
(2) So now I had to find the norms of large submatrices of A. Whilst I
remembered that the norm of a matrix is the largest amount by which it
stretches an element, I confess that for a short while I confused this idea
with that of the largest eigenvector of the matrix, and even began to
program with this in mind. The early results of the program alerted me to
the distinction between stretches that preserve direction and those that
don't. I looked up the result I needed, that the norm of a matrix M is the
square root of the largest eigenvalue of (MT)(M), where MT is the
transpose of M.
(3) To find the largest eigenvalue of a big matrix, I used the idea that if
a vector x is a linear combination of eigenvectors, say x = e_1 + e_2 + ...
+ e_n, then {[(MT)(M)]^k}(x) will be (lambda_1)^k.(e_1) +
(lambda_2)^k.(e_2) + ... + (lambda_n)^k.(e_n), so that the e_i with the
biggest lambda_i will increasingly dominate the vector formed by repeatedly
premultiplying x by (MT)(M). The UBASIC program below does exactly that.
After each iteration, it re-scales the vector to unit length by dividing it
by its length.
(4) The program generates a starting vector at random. I ran the program
with several values for the random seed to guard against unluckily getting
a starting vector expressible as a linear combination of the eigenvectors
_other_ than the one I wanted. Once the gap between successive estimates of
the largest eignevalue falls below the value set by "err", the program
doubles the size of the sub-matrix and starts again, using the estimated
eigenvector from the smaller matrix as the initial trial vector (early runs
of the program did not use this feature, as I wanted to check that there
were no big discontinuities in behaviour as the size of the matrix grew,
but for the big matirces needed to establish the 10 figures of accuracy, I
needed the time saving achieved by having a good first approximation).
(5) The answer!
In the output given in (7) below the gaps between norms of successive
matrices decrease by a factor of about 8 as matrix size doubles. Assuming
this pattern coninues, we can be sure of the answer to 10 figures: it
truncates to 1.274224152 and rounds to 1.274224153
(6) Here is the program:
1 'siam challenge problem 3
2 '
10 point 10 ' sets number of decimal places
20 Err=1/10^16 ' "error" between successive estimates
21 ' must become less than err before moving on to next matrix
30 Oldnorm=100
40 randomize 1
48 '
49 ' lines 60 and 70 set up data files used for storing vectors
50 ' File1 is the image after MT, file2 intermediate image after M.
55 ' The word values and the filling with random numbers in lines
56 ' 60 to 80 are not logical but UBASIC was giving me error
57 ' messages that didn't make sense either and this tinkering
58 ' seemed to cure the problem
59 '
60 open "si3vec1" as file1(7000) word 15
70 open "si3vec2" as file2(7000) word 15
80 for I=1 to 7000:file1(I)=rnd:file2(I)=rnd:next I
81 '
82 ' 90 and 100 set up files for saving output
83 '
90 open "si3res.txt" for output as #1
100 open "si3scre.txt" for output as #2
101 '
110 Siz=200
120 print #2,"screen output when initial size is";Siz
130 '
140 'fill vector1 at random
141 '
150 randomize
160 S=0
170 for I=1 to Siz
180 file1(I)=rnd
190 S+=abs(file1(I)^2)
200 next I
211 '
212 ' "Norm is a bit of a misnomer in 220. It is the norm of a vector
213 ' but on later passes through this code it's the square of the estimate
214 ' for the norm of the matrix.
216 '
220 Norm=sqrt(S)
249 '
250 'check for convergence
251 '
260 if abs(Norm-Oldnorm)<Err then Siz*=2:print "size=";Siz:print #2,"size=";Siz:print #1,Siz/2,sqrt(Norm):for I=Siz\2+1 to Siz:file1(I)=0:next I
270 Oldnorm=Norm
277 '
278 'rescale current vector to unit length
279 '
280 for I=1 to Siz
290 Z=file1(I)/Norm
300 file1(I)=Z
310 next I
311 '
317 ' The next two lines print garbage 1st time through but on later
318 ' passes print out the current estimate of the norm
319 '
320 print sqrt(Norm)
330 print #2,sqrt(Norm)
335 '
340 ' Multiply file1 By M and Store Result In file2
345 '
350 S=0
360 for I=1 to Siz
370 Newelt=0
380 for J=1 to Siz
390 Newelt+=(fnA(I,J)*file1(J))
400 next J
410 S+=abs(Newelt)^2
420 file2(I)=Newelt
430 next I
431 '
435 ' print the "stretch" from x to M(x) - another estimate of the norm
439 '
440 print "(";sqrt(S);")"
450 print #2"(";sqrt(S);")"
455 '
460 ' Multiply file2 By MT and Store Result In file1
465 '
470 S=0
480 for I=1 to Siz
490 Newelt=0
500 for J=1 to Siz
510 Newelt+=fnA(J,I)*file2(J)
520 next J
530 S+=abs(Newelt)^2
540 file1(I)=Newelt
550 next I
560 goto 250
797 '
798 'Routine for Calculating Elements Of M
799 '
800 fnA(I,J)
810 local N
820 N=I+J-2
830 return(1/(0.5*N*(N+1)+I))
(7) sample output:
screen output when initial size is 200
2.777128150913789805230078616851372410556616621222
( 0.125090649973425233649683956136349493058063328264 )
0.386069863653173486981635516867942626799516009373
( 1.273193932981053739736233891712848931482825390414 )
1.273702251258886387698802153871941570302882250772
( 1.274223441010932231642904044971638909178200320253 )
1.274223520179869730370254622280176495408698273634
( 1.274223601328154886373162525794155785813114414747 )
1.274223601340526627167487730515401658529319773736
( 1.274223601353207682816114789395583874821432626922 )
1.274223601353209616167334235143084955665039929703
( 1.274223601353211597855668952797578584886216043368 )
1.274223601353211598157796767105780367134016456156
( 1.274223601353211598467478297904121046936151927555 )
size= 400
1.274223601353211598467525511889486851556689268052
( 1.274223858394176413727028339127882946136219197314 )
1.274223969892378332360689005106707995586524521697
( 1.274224081392223999386371415817929529181804839855 )
1.274224081392226688701223845551857729333411253355
( 1.274224081392229441004820443587195221176362515875 )
1.274224081392229441397495365960643541330184262698
( 1.274224081392229441799988837978473663624745449699 )
size= 800
1.274224081392229441800050221819212161939090764464
( 1.274224114866918164430727075379474974886096625442 )
1.274224129284003428801755066287906259160099955518
( 1.274224143701141589768204819746910125030848979168 )
1.274224143701141641972818679651654423083391289476
( 1.274224143701141695375200070568882552200993031407 )
1.274224143701141695382661116097466311642403481075
( 1.274224143701141695390308719118469712673700073204 )
size= 1600
1.27422414370114169539030988548051231863663935352
( 1.274224147987001053758023288819531560887377992676 )
1.274224149826982798521302008681638616001916227161
( 1.27422415166696621795730058457566353577858707572 )
1.274224151666966218899601553222978298055515317455
( 1.274224151666966219863118862883905845291823395172 )
size= 3200
1.274224151666966219863250930801970451979602797302
( 1.2742241522101484486584930734180590928109778603 )
1.274224152443028887875950968641277888300854084173
( 1.274224152675909379710924267072870365374652416087 )
1.274224152675909379727105882597475168081641934905
( 1.274224152675909379743646033524044382302795308088 )
size= 6400
1.274224152675909379743648264126147888455923501466
( 1.274224152744342178641519689434302775887730671414 )
1.274224152773665789038775287233087330270597651217
( 1.274224152802989401083743928294146221069842644932 )
1.274224152802989401084012602383214855725106495874
( 1.27422415280298940108428715330417670103738065913 )
size= 12800
[The program elegantly crashes here as the size of the arrays is exceeded]