2 views (last 30 days)

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

David Goodmanson
on 23 May 2020

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.

David Goodmanson
on 23 May 2020

Hi ektor,

I was assuming that a_1 times h is the inner product (scalar product). that is, one number for the output, is that correct? Or do you mean element-by-element multiplicaton with no sum, so that you get an output of length(a_1) numbers?

David Goodmanson
on 3 Jun 2020

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.

per isakson
on 23 May 2020

Edited: per isakson
on 23 May 2020

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

per isakson
on 23 May 2020

Do you ask for why or how? I answer both.

"[...] constructing a complicated matrix A. The matrix A is a diagonal matrix" but I'm still not sure exactly what OP asked for.

Internally, Matlab stores arrays as linear sequences of elements together with some metadata, including the size of the array. Matlab uses column-major (and C row-major). Accessing array elements can be done with subscripts or linear indexing. (The functions, sub2ind and ind2sub convert between the two.) 1:N+1:end is the linear indicies of the diagonal elements.

David Goodmanson
on 24 May 2020

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 2 Comments

## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/530768-constructing-a-difficult-large-matrix#comment_858503

⋮## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/530768-constructing-a-difficult-large-matrix#comment_858503

## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/530768-constructing-a-difficult-large-matrix#comment_858618

⋮## Direct link to this comment

https://fr.mathworks.com/matlabcentral/answers/530768-constructing-a-difficult-large-matrix#comment_858618

Sign in to comment.