Suppose I have an NxN matrix A, an index vector V consisting of a subset of the numbers 1:N, and a value K, and I want to do this:
for i = V A(i,i) = K end Is there a way to do this in one statement w/ vectorization?
e.g. A(something) = K
The statement A(V,V) = K will not work, it assigns off-diagonal elements, and this is not what I want. e.g.:
>> A = zeros(5); >> V = [1 3 4]; >> A(V,V) = 1 A = 1 0 1 1 0 0 0 0 0 0 1 0 1 1 0 1 0 1 1 0 0 0 0 0 0 I usually use EYE for that:
A = magic(4) A(logical(eye(size(A)))) = 99 A = 99 2 3 13 5 99 10 8 9 7 99 12 4 14 15 99 Alternatively, you can just create the list of linear indices, since from one diagonal element to the next, it takes nRows+1 steps:
[nRows,nCols] = size(A); A(1:(nRows+1):nRows*nCols) = 101 A = 101 2 3 13 5 101 10 8 9 7 101 12 4 14 15 101 If you only want to access a subset of diagonal elements, you need to create a list of diagonal indices:
subsetIdx = [1 3]; diagonalIdx = (subsetIdx-1) * (nRows + 1) + 1; A(diagonalIdx) = 203 A = 203 2 3 13 5 101 10 8 9 7 203 12 4 14 15 101 Alternatively, you can create a logical index array using diag (works only for square arrays)
diagonalIdx = false(nRows,1); diagonalIdx(subsetIdx) = true; A(diag(diagonalIdx)) = -1 A = -1 2 3 13 5 101 10 8 9 7 -1 12 4 14 15 101 
diag first, before I remember to use eye
A(rowIdx, colIdx) = -1 assigns the value to all combinations of rowIdx/colIdx. Thus, if you want to assign multiple elements of an array in one go, linear indices are the way to go.sub2ind which returns you the linear index for given rowIdx/colIdxsub2ind would indeed make sense that way.>> tt = zeros(5,5) tt = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> tt(1:6:end) = 3 tt = 3 0 0 0 0 0 3 0 0 0 0 0 3 0 0 0 0 0 3 0 0 0 0 0 3 and more general:
>> V=[1 2 5]; N=5; >> tt = zeros(N,N); >> tt((N+1)*(V-1)+1) = 3 tt = 3 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 This is based on the fact that matrices can be accessed as one-dimensional arrays (vectors), where the 2 indices (m,n) are replaced by a linear mapping m*N+n.
>> B=[0,4,4;4,0,4;4,4,0] B = 0 4 4 4 0 4 4 4 0 >> v=[1,2,3] v = 1 2 3 >> B(eye(size(B))==1)=v %insert values from v to eye positions in B B = 1 4 4 4 2 4 4 4 3 
A = zeros(7,6); V = [1 3 5]; [n m] = size(A); diagIdx = 1:n+1:n*m; A( diagIdx(V) ) = 1 A = 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 I'd use sub2ind and pass the diagonal indices as both x and y parameters:
A = zeros(4) V=[2 4] idx = sub2ind(size(A), V,V) % idx = [6, 16] A(idx) = 1 % A = % 0 0 0 0 % 0 1 0 0 % 0 0 0 0 % 0 0 0 1 Suppose K is the value. The command
A=A-diag(K-diag(A)) may be a bit faster
>> A=randn(10000,10000); >> tic;A(logical(eye(size(A))))=12;toc Elapsed time is 0.517575 seconds.
>> tic;A=A+diag((99-diag(A)));toc Elapsed time is 0.353408 seconds.
But it consumes more memory.
A(logical(eye(size(A))))=K flexible fast and reliableI use this small inline function in finite difference code.
A=zeros(6,3); range=@(A,i)[1-min(i,0):size(A,1)-max(i+size(A,1)-size(A,2),0 ) ]; Diag=@(A,i) sub2ind(size(A), range(A,i),range(A,i)+i ); A(Diag(A, 0))= 10; %set diagonal A(Diag(A, 1))= 20; %equivelent to diag(A,1)=20; A(Diag(A,-1))=-20; %equivelent to diag(A,-1)=-20; It can be easily modified to work on a sub-range of the diagonal by changing the function range.
