sparse matrix 3 diagonals and 3 elements
wma
i need to create a sparse matrix with 3 diagonals. Each diagonal stores only one element.
My matrix looks like this:
| a 0 . 0 b 0 . 0 c 0 . 0|
| 0 a . 0 0 b . 0 0 c . 0|
| . . . 0 . . . . 0 . . . 0|
| 0 0 0 a 0 0 0 b 0 0 0 c|
i found a function called spdiags and stuck with it.
can anyone help please?
Jos (10584)
You can simple concatenate the individual diagonal sparse matrices. Here is an example:
A = 10 ; B = 20 ;
Siz = 4 ;
s1 = spdiags(repmat(A,Siz,1),0,Siz,Siz)
s2 = spdiags(repmat(B,Siz,1),0,Siz,Siz)
result = [s1 s2]
full(result) % nicer display
hth
Jos
wma
>
> A = 10 ; B = 20 ;
> Siz = 4 ;
>
> s1 = spdiags(repmat(A,Siz,1),0,Siz,Siz)
> s2 = spdiags(repmat(B,Siz,1),0,Siz,Siz)
> result = [s1 s2]
> full(result) % nicer display
>
> hth
> Jos
Thank you for the big help!
John D'Errico
So what are you stuck about? spdiags is a very
reasonable and proper solution here.
It appears that the matrix you wish to create
is basically three identity matrices, scaled and
joined together by horizontal concatenation.
Assuming that to be an accurate description,
how would you write the spdiags call? TRY IT!
Try it yourself before you read what I will write
below, else you will not learn as much.
...
...
...
...
...
...
Did you try it? Try again. Yes, I'm being picky
here. But you will learn only by looking at the
help, and making an attempt of your own.
...
...
...
...
...
...
I give up. Lets see how you might build this
matrix.
a = 1;
b = 2;
c = 3;
% assume nxn blocks here.
n = 100000;
% a simple solution, using brute force
tic
M0 = [speye(n)*a,speye(n)*b,speye(n)*c];
toc
% Elapsed time is 0.287485 seconds.
% A call to spy shows that this does as expected
spy(M0)
% Can we use spdiags here though?
tic
M1 = spdiags(repmat([a,b,c],n,1),[0,n,2*n],n,3*n);
toc
% Elapsed time is 0.564054 seconds.
Surprisingly, the single call to spdiags took longer.
Yes, I was surprised too. Can we do any better than
the speye calls? Perhaps.
tic
M2 = sparse(repmat((1:n)',1,3), ...
bsxfun(@plus,(1:n)',n*(0:2)), ...
repmat([a,b,c],n,1),n,3*n);
toc
% Elapsed time is 0.256570 seconds.
Are these matrices all the same? Yes.
nnz(M0-M2)
ans =
0
nnz(M0-M1)
ans =
0
HTH,
John
wma
M0 = [speye(n)*a,speye(n)*b,speye(n)*c];
this solution is fastest for me.
This variant:
M2 = sparse(repmat((1:n)',1,3), ...
bsxfun(@plus,(1:n)',n*(0:2)), ...
repmat([a,b,c],n,1),n,3*n);
is suddently not as fast as you posted.
How does matlab handle the internal data?
Do i have an Array of 3*n elements or only 3 elements as the most optimal solution?
Is there a possibility to get the Matrix size in bytes like C++ sizeof(X) ?
John D'Errico
> Thanks for the idea with:
>
> M0 = [speye(n)*a,speye(n)*b,speye(n)*c];
>
> this solution is fastest for me.
>
> This variant:
> M2 = sparse(repmat((1:n)',1,3), ...
> bsxfun(@plus,(1:n)',n*(0:2)), ...
> repmat([a,b,c],n,1),n,3*n);
> is suddently not as fast as you posted.
Different methods may vary in speed, depending on
the computer and matlab release you are running.
> How does matlab handle the internal data?
> Do i have an Array of 3*n elements or only 3 elements as the most optimal solution?
How can/should matlab know that these elements
are exactly equal? It would take far more effort
and time to determine this fact than the minor
amount of memory saved is worth. Given that
your specific case is a rare one, there is no need
to optimize for such an uncommon case.
> Is there a possibility to get the Matrix size in bytes like C++ sizeof(X) ?
help whos
John
wma
> are exactly equal? It would take far more effort
> and time to determine this fact than the minor
> amount of memory saved is worth. Given that
> your specific case is a rare one, there is no need
> to optimize for such an uncommon case.
Thank you for your time to answer me but at this point you are maybe unright. I am writing my theisis about TLS problems in Gauss-Helmert Modell. It is the most common case in TLS adjustment problems. In most of B(Design)-Matrices we have such a case.
| a 0 0 b 0 0 c 0 0 |
| 0 a 0 0 b 0 0 c 0 | = B
| 0 0 a 0 0 b 0 0 c |
such matrices are growing with the count of points and if we have data from terrestrial laserscanners, some point clouds can reach 30M points and more,
in which case B*B' and inv(B*B') are unfortunally big and have only one value.
I know many situations outside of TLS adjustment where matlab can benefit from singlevalued arrays.
It is not needed to check this cases automatically, its enough to have such functions.
John D'Errico
YOU would find this to be of value. But the fact is,
if would give little gain to most of the rest of the
world, and not enough gain to justify a highly
specific matrix class that will be difficult to specify.
The fact is you CAN do almost anything you need
by simple linear algebra. You define the array B as
the concatenation of three identity matrices, each
times a different scalar. But any matrix multiply
can then be broken down using blocked matrices.
So you can indeed do what you wish already.
Finally since you do feel to need this capability in
MATLAB, you can build it yourself, quite trivially so
as a new class in MATLAB.
So where is the problem?
John