{"id":8810,"date":"2022-07-20T23:38:23","date_gmt":"2022-07-21T03:38:23","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=8810"},"modified":"2022-07-20T23:45:59","modified_gmt":"2022-07-21T03:45:59","slug":"an-interesting-and-perhaps-new-matrix","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2022\/07\/20\/an-interesting-and-perhaps-new-matrix\/","title":{"rendered":"An Interesting, and Perhaps New, Matrix"},"content":{"rendered":"\r\n\r\n<div class=\"content\"><!--introduction--><p>I do not recall having seen this matrix before, but I will not be surprised to learn that somebody else already knows all about it, especially if that person's name is Nick.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#17afe407-67a2-462b-a10b-df348fccb837\"><tt>Q<\/tt><\/a><\/li><li><a href=\"#044996c2-78b8-4a23-a809-6ed8ab3e3607\"><tt>D<\/tt><\/a><\/li><li><a href=\"#eca1bb02-1ba1-4bbf-9ed0-39c244e018eb\">O.E.I.S.<\/a><\/li><li><a href=\"#b7778074-8d3f-46ad-858b-13f27d2f99bd\"><tt>R<\/tt><\/a><\/li><li><a href=\"#1e78f08f-c39d-4810-93bb-b0428456710c\">Condition<\/a><\/li><li><a href=\"#c60c371e-4378-4499-8242-376b1b37da8c\">Extra Credit<\/a><\/li><\/ul><\/div><h4><tt>Q<\/tt><a name=\"17afe407-67a2-462b-a10b-df348fccb837\"><\/a><\/h4><p>I've been investigating the matrices generated by this elegant one-liner.<\/p><pre class=\"codeinput\">    Q = @(n) (-n:n).^2 + (-n:n)'.^2;\r\n<\/pre><p>The <tt>Q<\/tt> is for \"quadratic\".<\/p><p>The middle column contains the squares of the integers from <tt>-n<\/tt> to <tt>n<\/tt>. So does the middle row. The apostrophe summons singleton expansion. The resulting matrix has order <tt>2*n+1<\/tt>. Here is <tt>Q(5)<\/tt>.<\/p><pre class=\"codeinput\">    Q5 = Q(5)\r\n<\/pre><pre class=\"codeoutput\">\r\nQ5 =\r\n\r\n    50    41    34    29    26    25    26    29    34    41    50\r\n    41    32    25    20    17    16    17    20    25    32    41\r\n    34    25    18    13    10     9    10    13    18    25    34\r\n    29    20    13     8     5     4     5     8    13    20    29\r\n    26    17    10     5     2     1     2     5    10    17    26\r\n    25    16     9     4     1     0     1     4     9    16    25\r\n    26    17    10     5     2     1     2     5    10    17    26\r\n    29    20    13     8     5     4     5     8    13    20    29\r\n    34    25    18    13    10     9    10    13    18    25    34\r\n    41    32    25    20    17    16    17    20    25    32    41\r\n    50    41    34    29    26    25    26    29    34    41    50\r\n\r\n<\/pre><p>I like the contour plot.<\/p><pre class=\"codeinput\">    contourf(Q(100))\r\n    axis <span class=\"string\">square<\/span>\r\n    colorbar\r\n    title(<span class=\"string\">'Q(100)'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/disc_blog_01.png\" alt=\"\"> <h4><tt>D<\/tt><a name=\"044996c2-78b8-4a23-a809-6ed8ab3e3607\"><\/a><\/h4><p>For another blog post under development, I need a logical mask that carves a circular region out of graphic.  This disc does the job.<\/p><pre class=\"codeinput\">    D = @(n) Q(n) &lt;= n^2;\r\n<\/pre><p>Here is my carver.<\/p><pre class=\"codeinput\">    spy(D(100))\r\n    title(<span class=\"string\">'D(100)'<\/span>)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/disc_blog_02.png\" alt=\"\"> <p>Did you notice the digits in the count of nonzeros in <tt>D(100)<\/tt>? It happens whenever <tt>n<\/tt> is a power of 10.<\/p><pre class=\"codeinput\">    fprintf(<span class=\"string\">'%15s %12s\\n'<\/span>,<span class=\"string\">'n'<\/span>,<span class=\"string\">'nnz(D(n))'<\/span>)\r\n    <span class=\"keyword\">for<\/span> n = 10.^(0:4)\r\n        fprintf(<span class=\"string\">'%15d %12d\\n'<\/span>,n, nnz(D(n)))\r\n    <span class=\"keyword\">end<\/span>\r\n<\/pre><pre class=\"codeoutput\">              n    nnz(D(n))\r\n              1            5\r\n             10          317\r\n            100        31417\r\n           1000      3141549\r\n          10000    314159053\r\n<\/pre><h4>O.E.I.S.<a name=\"eca1bb02-1ba1-4bbf-9ed0-39c244e018eb\"><\/a><\/h4><p>A classic problem, described in the <a href=\"https:\/\/oeis.org\/A000328\">Online Encyclopedia of Integer Sequences<\/a>, asks how many points with integer coordinates lie within a disc of increasing radius. Our nonzero count provides the answer.<\/p><pre class=\"codeinput\">    fprintf(<span class=\"string\">'%15s %8s\\n'<\/span>,<span class=\"string\">'n'<\/span>,<span class=\"string\">'a(n)'<\/span>)\r\n    <span class=\"keyword\">for<\/span> n = [1:15  99:101  499:501  999:1001]\r\n        <span class=\"keyword\">if<\/span> mod(n,100) == 99\r\n            fprintf(<span class=\"string\">'%15s %8s\\n'<\/span>,<span class=\"string\">'-'<\/span>,<span class=\"string\">'-'<\/span>)\r\n        <span class=\"keyword\">end<\/span>\r\n        a(n) = nnz(D(n));\r\n        fprintf(<span class=\"string\">'%15d %8d\\n'<\/span>,n,a(n))\r\n    <span class=\"keyword\">end<\/span>\r\n<\/pre><pre class=\"codeoutput\">              n     a(n)\r\n              1        5\r\n              2       13\r\n              3       29\r\n              4       49\r\n              5       81\r\n              6      113\r\n              7      149\r\n              8      197\r\n              9      253\r\n             10      317\r\n             11      377\r\n             12      441\r\n             13      529\r\n             14      613\r\n             15      709\r\n              -        -\r\n             99    30757\r\n            100    31417\r\n            101    32017\r\n              -        -\r\n            499   782197\r\n            500   785349\r\n            501   788509\r\n              -        -\r\n            999  3135157\r\n           1000  3141549\r\n           1001  3147833\r\n<\/pre><h4><tt>R<\/tt><a name=\"b7778074-8d3f-46ad-858b-13f27d2f99bd\"><\/a><\/h4><p>Taking the reciprocals of the matrix entries, and reducing the range of the anonymous index, produces a matrix that behaves a bit like the Hilbert matrix, <tt>hilb(n)<\/tt>.<\/p><pre class=\"codeinput\">    R = @(n) 1.\/((1:n).^2 + (1:n)'.^2);\r\n<\/pre><p>Here are the 5-by-5's.<\/p><pre class=\"codeinput\">    format <span class=\"string\">rat<\/span>\r\n    R5 = R(5)\r\n    H5 = hilb(5)\r\n<\/pre><pre class=\"codeoutput\">\r\nR5 =\r\n\r\n       1\/2            1\/5            1\/10           1\/17           1\/26    \r\n       1\/5            1\/8            1\/13           1\/20           1\/29    \r\n       1\/10           1\/13           1\/18           1\/25           1\/34    \r\n       1\/17           1\/20           1\/25           1\/32           1\/41    \r\n       1\/26           1\/29           1\/34           1\/41           1\/50    \r\n\r\n\r\nH5 =\r\n\r\n       1              1\/2            1\/3            1\/4            1\/5     \r\n       1\/2            1\/3            1\/4            1\/5            1\/6     \r\n       1\/3            1\/4            1\/5            1\/6            1\/7     \r\n       1\/4            1\/5            1\/6            1\/7            1\/8     \r\n       1\/5            1\/6            1\/7            1\/8            1\/9     \r\n\r\n<\/pre><h4>Condition<a name=\"1e78f08f-c39d-4810-93bb-b0428456710c\"><\/a><\/h4><p>Going away from the diagonal, the elements of <tt>R(n)<\/tt> decay more rapidly than those of <tt>hilb(n)<\/tt>, so <tt>R(n)<\/tt> is better conditioned than <tt>hilb(n)<\/tt>.<\/p><pre class=\"codeinput\">    format <span class=\"string\">short<\/span> <span class=\"string\">e<\/span>\r\n    fprintf(<span class=\"string\">'%15s %12s %12s\\n'<\/span>,<span class=\"string\">'n'<\/span>,<span class=\"string\">'cond R'<\/span>,<span class=\"string\">'cond H'<\/span>)\r\n    <span class=\"keyword\">for<\/span> n = 1:12\r\n        fprintf(<span class=\"string\">'%15d %12.2e %12.2e\\n'<\/span>,n,cond(R(n)),cond(hilb(n)))\r\n    <span class=\"keyword\">end<\/span>\r\n<\/pre><pre class=\"codeoutput\">              n       cond R       cond H\r\n              1     1.00e+00     1.00e+00\r\n              2     1.53e+01     1.93e+01\r\n              3     2.04e+02     5.24e+02\r\n              4     2.59e+03     1.55e+04\r\n              5     3.21e+04     4.77e+05\r\n              6     3.89e+05     1.50e+07\r\n              7     4.67e+06     4.75e+08\r\n              8     5.54e+07     1.53e+10\r\n              9     6.53e+08     4.93e+11\r\n             10     7.65e+09     1.60e+13\r\n             11     8.92e+10     5.22e+14\r\n             12     1.04e+12     1.62e+16\r\n<\/pre><h4>Extra Credit<a name=\"c60c371e-4378-4499-8242-376b1b37da8c\"><\/a><\/h4><p>What is the rank of <tt>Q(n)<\/tt>?  Why?  See <a href=\"https:\/\/epubs.siam.org\/doi\/epdf\/10.1137\/20M1358694\">a paper in SIAM Review<\/a> by Strang and Moler.<\/p><p>Why is the table of values for <tt>nnz(D(10^k))<\/tt> so short? How might you extend this table?<\/p><p>Investigate <tt>R(n)<\/tt>.  Is it positive definite?  What are its eigenvalues?  What is its inverse?  What is the sign pattern of the elements of its inverse?  For what values of <tt>n<\/tt> can you compute the inverse reliably using floating point arithmetic? How does all this compare with <tt>hilb(n)<\/tt> and <tt>invhilb(n)<\/tt> ?<\/p><p>Comments in the Comments, or in email to me, are welcome.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_4c513405f3644e5fbcae1bf70683cb1d() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='4c513405f3644e5fbcae1bf70683cb1d ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 4c513405f3644e5fbcae1bf70683cb1d';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2022 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_4c513405f3644e5fbcae1bf70683cb1d()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2022a<br><\/p><\/div><!--\r\n4c513405f3644e5fbcae1bf70683cb1d ##### SOURCE BEGIN #####\r\n%% An Interesting, and Perhaps New, Matrix\r\n% I do not recall having seen this matrix before, but\r\n% I will not be surprised to learn that somebody else already knows \r\n% all about it, especially if that person's name is Nick.\r\n\r\n%% |Q| \r\n%\r\n% I've been investigating the matrices generated by this elegant one-liner.\r\n \r\n    Q = @(n) (-n:n).^2 + (-n:n)'.^2;\r\n\r\n%%\r\n% The |Q| is for \"quadratic\". \r\n%\r\n% The middle column contains the squares of the integers from |-n| to |n|.\r\n% So does the middle row. The apostrophe summons singleton expansion.\r\n% The resulting matrix has order |2*n+1|.\r\n% Here is |Q(5)|.\r\n\r\n    Q5 = Q(5)\r\n\r\n%%\r\n% I like the contour plot.\r\n\r\n    contourf(Q(100))\r\n    axis square\r\n    colorbar\r\n    title('Q(100)')\r\n\r\n%% |D|\r\n% For another blog post under development, I need a logical mask that\r\n% carves a circular region out of graphic.  This disc does the job.\r\n \r\n    D = @(n) Q(n) <= n^2;\r\n\r\n%%\r\n% Here is my carver.\r\n\r\n    spy(D(100))\r\n    title('D(100)')\r\n\r\n%% \r\n% Did you notice the digits in the count of\r\n% nonzeros in |D(100)|?\r\n% It happens whenever |n| is a power of 10.\r\n\r\n    fprintf('%15s %12s\\n','n','nnz(D(n))')\r\n    for n = 10.^(0:4)\r\n        fprintf('%15d %12d\\n',n, nnz(D(n)))\r\n    end\r\n\r\n%% O.E.I.S.\r\n% A classic problem, described in the <https:\/\/oeis.org\/A000328\r\n% Online Encyclopedia of Integer Sequences>, asks how many points with\r\n% integer coordinates lie within a disc of increasing radius.\r\n% Our nonzero count provides the answer.\r\n\r\n    fprintf('%15s %8s\\n','n','a(n)')\r\n    for n = [1:15  99:101  499:501  999:1001]\r\n        if mod(n,100) == 99\r\n            fprintf('%15s %8s\\n','-','-')\r\n        end\r\n        a(n) = nnz(D(n));\r\n        fprintf('%15d %8d\\n',n,a(n))\r\n    end  \r\n\r\n%% |R|\r\n% Taking the reciprocals of the matrix entries, and reducing the range of\r\n% the anonymous index, produces a matrix that behaves a bit like\r\n% the Hilbert matrix, |hilb(n)|.\r\n\r\n    R = @(n) 1.\/((1:n).^2 + (1:n)'.^2);\r\n\r\n%%\r\n% Here are the 5-by-5's.\r\n\r\n    format rat\r\n    R5 = R(5)\r\n    H5 = hilb(5)\r\n\r\n%% Condition\r\n% Going away from the diagonal,\r\n% the elements of |R(n)| decay more rapidly than those of |hilb(n)|,\r\n% so |R(n)| is better conditioned than |hilb(n)|.\r\n\r\n    format short e\r\n    fprintf('%15s %12s %12s\\n','n','cond R','cond H')\r\n    for n = 1:12\r\n        fprintf('%15d %12.2e %12.2e\\n',n,cond(R(n)),cond(hilb(n)))\r\n    end\r\n\r\n%% Extra Credit\r\n% What is the rank of |Q(n)|?  Why?  See \r\n% <https:\/\/epubs.siam.org\/doi\/epdf\/10.1137\/20M1358694\r\n% a paper in SIAM Review> by Strang and Moler.\r\n%\r\n% Why is the table of values for |nnz(D(10^k))| so short?\r\n% How might you extend this table?\r\n%\r\n% Investigate |R(n)|.  Is it positive definite?  What are its\r\n% eigenvalues?  What is its inverse?  What is the sign pattern of\r\n% the elements of its inverse?  For what values of |n| can you\r\n% compute the inverse reliably using floating point arithmetic?\r\n% How does all this compare with |hilb(n)| and |invhilb(n)| ?\r\n%\r\n% Comments in the Comments, or in email to me, are welcome. \r\n\r\n##### SOURCE END ##### 4c513405f3644e5fbcae1bf70683cb1d\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/disc_blog_00.png\" class=\"img-responsive attachment-post-thumbnail size-post-thumbnail wp-post-image\" alt=\"\" decoding=\"async\" loading=\"lazy\" \/><\/div><!--introduction--><p>I do not recall having seen this matrix before, but I will not be surprised to learn that somebody else already knows all about it, especially if that person's name is Nick.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2022\/07\/20\/an-interesting-and-perhaps-new-matrix\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":8828,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[5,4,6,16],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/8810"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=8810"}],"version-history":[{"count":3,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/8810\/revisions"}],"predecessor-version":[{"id":8834,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/8810\/revisions\/8834"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media\/8828"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=8810"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=8810"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=8810"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}