# Computing norm of a matrix with positive entries

Is there anything in Mathematica that can help me compute norm of a matrix with non-negative entries faster than Norm? Approach below works for matrix[1000], but too slow for matrix[100000]

matrix[d_] := Normalize /@ UpperTriangularize[ConstantArray[1., {d, d}]];
Norm@matrix[1000]


From the math here, I expect Norm@matrix[10^5] to be close to 262.99434, was curious to check how large approximation error is

predictedNorm[d_] := (
phi[\[Lambda]_] := Sum[(-\[Lambda])^n/(n!)^2, {n, 0, Infinity}];
l0 = \[Lambda] /.
First@FindRoot[phi[\[Lambda]] == 0, {\[Lambda], 1.4},
WorkingPrecision -> 50, AccuracyGoal -> 20,
PrecisionGoal -> 20];
Sqrt[d/l0]
);
predictedNorm[100000] (* 262.99434901 *)

• Since your triangular matrix M is dense enough, it is already very memory-demanding to store 100000 x 100000 matrix elements of MM* which are all needed to compute the norm. Since you already know the kernel function (in Math forum), I advise you to construct "linear transformation" M and M* acting on 100000 dim vector space, which is of managable size. Then, the norm of M, which is the largest eigenvalue of MM*, is well estimated from the growth rate of repeated application of M and M* on random initial vector. Commented Jul 13 at 2:37
• (Now I see this...) Right, what @A.Kato said. Commented Jul 13 at 3:01

You're not going top fit that matrix just anywhere. You need a specialized matrix-times-vector routine (and vice versa, for the case at hand). I won't show the simple confirmation tests I did, but these do the job.

matTimesVec[vec_] := Module[{dim = Length[vec]},
Table[Total[vec[[j ;;]]]/Sqrt[(dim + 1. - j)]
, {j, dim}]]

vecTimesMat[vec_] := Module[{dim = Length[vec]},
Table[vec[[1 ;; j]] . (1/Sqrt[(dim + 1. - Range[j])])
, {j, dim}]]


Now we do power iterations, in effect computing, for a random vector, updates of the form vec = Transpose[mat].mat.vec, using the fact that the transpose can instead be done as vec.mat. I've found this seems to stqabilize after a few iterations so I do 10 such. Then the norm is the square root of the penultimate vector over the last vector, in all positions not driven to zero.

In[159]:= siz = 10^4;
AbsoluteTiming[
rvec = RandomReal[1, siz];
Do[vlast = rvec;
rvec = matTimesVec[rvec];
rvec = vecTimesMat[rvec],
10];
Sqrt[rvec[[1 ;; 5]]/vlast[[1 ;; 5]]]
]

(* Out[160]= {3.9033, {83.1682, 83.1682, 83.1682, 83.1682, 83.1682}} *)


So 4 seconds to handle 10000. It's an O(n^2) algorithm so give several minutes for 10^5.

siz = 10^5;
Timing[
rvec = RandomReal[1, siz];
Do[vlast = rvec;
rvec = matTimesVec[rvec];
rvec = vecTimesMat[rvec],
10];
Sqrt[rvec[[1 ;; 5]]/vlast[[1 ;; 5]]]
]

(* {356.045, {262.995, 262.995, 262.995, 262.995, 262.995}} *)