constructing a difficult large matrix

Question

ektor 2020-5-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/530768-constructing-a-difficult-large-matrix

评论： David Goodmanson 2020-6-3

采纳的回答： David Goodmanson

在 MATLAB Online 中打开

Dear all,

I am interested in constructing a complicated matrix A.

The matrix A is a diagonal matrix, where each diagonal element 'a_i' is of 1X1000, having the following patter

a_1=[p 1 0 0 0 0 0 .....0 ];
a_2=[p^2 p 1 0 0 0..... 0];
a_2=[p^3 p^2 p 1 0 0......0];
a_2=[p^4 p^3 p^2 p 1 0 0......0 ];
a_2=[p^5 p^4 p^3 p^2 p 1 0 ...0];
.
.
.
a_1000=[p^1000 p^999 p^998.........1];

And I want to multiply each element "a_i" with the vector h, where h is of dimension 1000X1, which contains just numbers of no particular pattern

Is there any way of doing that fast?

Maybe by using sparse or speye?

Thank you

2 个评论
显示无隐藏无

per isakson 2020-5-23

"multiply each element "a_i" with the vector h, where h is of dimension 1000X1" the product will that be a scalar, vector or matrix?

"doing that fast" how fast is "fast" ?

ektor 2020-5-23

Hi isakson,

It will be a scalar. So each element of a_i will be multiplied by the corresponding element of h to produce a scalar (like inner product). Sorry for not clarifying first.

When I say fast, I mean not time consuming.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Goodmanson 2020-5-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/530768-constructing-a-difficult-large-matrix#answer_436833

在 MATLAB Online 中打开

Hi ektor,

I presume you are interested not so much in the in the matrix A as in the 1000 scalars you get after you multiply each row of A by the vector h. To make things a little easier I included an extra row a_0 at the beginning.

a_0=[1 0 0 0 0 0 0 .....0 ];
a_1=[p 1 0 0 0 0 0 .....0 ];
a_2=[p^2 p 1 0 0 0..... 0];
a_3=[p^3 p^2 p 1 0 0......0];
a_4=[p^4 p^3 p^2 p 1 0 0......0 ];
a_5=[p^5 p^4 p^3 p^2 p 1 0 ...0];
.
.
a_1000=[p^1000 p^999 p^998.........1];

then

Ah = filter(1,[1 -p],h)

is the same as A*h. It's very fast. There is one exra element at the beginning corresponding to row a_0 but you can always delete it.

5 个评论
显示 3更早的评论隐藏 3更早的评论

ektor 2020-5-23

HI david,

Just a quick question: " There is one exra element at the beginning corresponding to row a_0 but you can always delete it. "

what does this mean?

David Goodmanson 2020-6-3

Hi ektor,

I meant that your matrix has 1000 rows and 1001 columns. The output is the inner product of h with each row, so 1000 values in all. However, for the filter function to reproduce that result, it was convenient to add a row at the top that I called a0. Then there are 1001 rows and 1001 values in the result. But the first result is due to the added row, so if you delete it you are back to the 1000 values in the original problem.

请先登录，再进行评论。

Answer 2

per isakson 2020-5-23

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/530768-constructing-a-difficult-large-matrix#answer_436883

编辑：per isakson 2020-5-23

在 MATLAB Online 中打开

Here are three functions, two of which are based on the answer of David Goodmanson. I think all of them are fast enough. However, more important than speed is that their results are correct. More test are needed.

A small test with numbers that I'm able to check.

>> cssm(8)
ans =
     3     0     0     0     0     0     0     0
     0     7     0     0     0     0     0     0
     0     0    15     0     0     0     0     0
     0     0     0    31     0     0     0     0
     0     0     0     0    63     0     0     0
     0     0     0     0     0   127     0     0
     0     0     0     0     0     0   255     0
     0     0     0     0     0     0     0   510
>> 
>> tic,cssm(1e4);toc
Elapsed time is 1.155732 seconds.
>> tic,M=dgsp(1e4);toc
Elapsed time is 0.490225 seconds.
>> tic,M=dg(1e4);toc
Elapsed time is 0.014082 seconds.
function    M = cssm( N )
    %%
    h = ones( N, 1 );
    p = 2;
    pr = 1;
    a = zeros( N, N );
    a( N+1 : (N+1) : end ) = 1;
    for jj = 1 : N
        pr = pr * p;
        a( jj : (N+1) : end-(jj-2)*N ) = pr;
    end
    
    M = zeros( N, N );
    for jj = 1 : N
        M(jj,jj) = a(jj,:) * h; 
    end
end
function    M = dgsp( N )
    %%
    h   = ones( N, 1 );
    tmp = [ 1; h ];
    p = 2;
    %%
    Ah = filter( 1, [1,-p], tmp );
    M = sparse( zeros( N, N ) );    % better: sparse(10,10,0); 
    M( 1 : N+1 : end ) = Ah( 2 : end );
end
function    M = dg( N )
    %%
    h   = ones( N, 1 );
    tmp = [ 1; h ];
    p = 2;
    %%
    Ah = filter( 1, [1,-p], tmp );
    M = zeros( N, N );
    M( 1 : N+1 : end ) = Ah( 2 : end );
end

4 个评论
显示 2更早的评论隐藏 2更早的评论

David Goodmanson 2020-5-24

I'm fine with the how, I was inquiring about the why, and now I see the direction you took. In the statement "A is a diagonal matrix, where each diagonal element 'a_i' is of 1X1000", if A is interpreted to be a very large block diagonal matrix, you do end up with the large diagonal matrix you have.