I have the below command in Matlab, I could not find how to do it in Julia. Any suggestion please?
sparse((1:5)', [11,22,33,44,55], [1,2,3,4,5])
1 Like
From help?> sparse
sparse(I, J, V,[ m, n, combine])
Create a sparse matrix S of dimensions m x n such that S[I[k], J[k]] = V[k].
Is this what you want?
julia> sparse(1:5, 1:5, [11,22,33,44,55])
5Γ5 SparseMatrixCSC{Int64, Int64} with 5 stored entries:
11 β
β
β
β
β
22 β
β
β
β
β
33 β
β
β
β
β
44 β
β
β
β
β
55
If so, maybe Diagonal would be better choice:
julia> Diagonal([11,22,33,44,55])
5Γ5 Diagonal{Int64, Vector{Int64}}:
11 β
β
β
β
β
22 β
β
β
β
β
33 β
β
β
β
β
44 β
β
β
β
β
55
2 Likes
In that case, itβs even closer to what you do in MATLAB, Iβve just abbreviated your [1,2,3,4,5] to 1:5, but it could be as an array like you have:
julia> sparse(1:5,[11,22,33,44,55],1:5)
5Γ55 SparseMatrixCSC{Int64, Int64} with 5 stored entries:
β β β β β β β β β β β β β β β β β β β β β β’β β β β β β
β β β β β β β β β β β β β β β β β β β β β β β β β β β β
julia> Array(ans)
5Γ55 Matrix{Int64}:
0 0 0 0 0 0 0 0 0 0 1 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 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 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 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 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 0 0 0 0 4 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5
2 Likes
@Jeff_Emanuel Thank you very much! Based on this, I have a related issue.
I create matrix F and then needs to be updated in a for loop as seen below. However, this way is time consuming (due the creation of sparseMatrix at each update). Is there another faster way?
Nh = [11,22,33,44,55]; # fixed
Ik = [10,12,23,34,45]; # variable
F = sparse(1:5,Nh,Ik)
iter=0;
for _ in 1:100
iter += 1;
Ik = ik*2*iter;
F = F.*iter + sparse(1:5,Nh,Ik)
end
Iβm getting the following error:
ERROR: LoadError: UndefVarError: ik not defined
Does ik change every iteration? Edit: ah maybe a type and it should be IK?
1 Like
Youβre always updating the same entries of your sparse matrix, so it would be vastly more efficient to update the entries in-place by mutating F.nzval.
Your code doesnβt make much sense because (even with the ik βΆ Ik typo fixed) it yields a zero matrix F the end because the integer entries overflow (perhaps you wanted a floating-point matrix instead? unlike Matlab, Julia numbers do not all default to floating-point values), but here is a similar example that updates the entries by adding random numbers:
irow = 1:5 # row indices
icol = [11,22,33,44,55]; # col indices
F = sparse(irow, icol, rand(5))
for iter = 1:100
@. F.nzval = iter * F.nzval + rand() # update nonzero entries in-place
end
which gives something like:
julia> @show F;
F = sparse([1, 2, 3, 4, 5], [11, 22, 33, 44, 55], [1.496910562914082e158, 1.8181168960906027e158, 8.993259798670023e157, 1.4533474201806958e158, 1.347203531062703e158], 5, 55)
3 Likes
I am sorry, yes it should be IK
It seems like you are only changing the data entries of the sparse matrix during the loop. Hence, you can precompute the data entries, IK, and then create the sparse matrix.
Simplifying your math in the code, we can the following simplification of the for loop:
Ik = [10.0, 12.0, 23.0, 34.0, 45.0] # variable
Ik .*= 2.0^100 * 2 * prod(1:100.0)
F = sparse(1:5, 11:11:55, Ik)
Timing gives the following:
129.232 ΞΌs (2447 allocations: 354.84 KiB) # Original code
560.737 ns (13 allocations: 1.80 KiB) # The simplified code above
Try to obtain the values of the final sparse matrix and then create a sparse matrix. And when obtaining the values of the final sparse matrix, simplify the mathematical expression.
Note: as your code is posted you will get integer overflow and the resulting matrix will be 0s. Using prod(1:100.0) instead of prod(1:100), 2.0^100 instead of 2^100, and Ik = [10.0, β¦] instead of [10, β¦] ensures that floats are used and thus the product can be computed successfully, without overflow. (Floats can represent larger numbers in memory compared to Ints)
Since Nh is fixed you can also use a UnitRange for it instead of allocating memory for the array. [11, 22, 33, 44, 55] turns into 11:11:55.
2 Likes
Thank you very much for your replies. I tried to simplify my real problem but it seems I have to mention it exactly as I am afraid I am missing important part. Below is exactly as from my code:
using SparseArrays, LinearAlgebra
R = sparse([zeros(1,8); Matrix(1.0I, 8-1, 8-1) zeros(8-1,1)]);
Nh = Int8[6, 6, 6, 8, 7, 7]; # fixed
Ik = [0, 1, -4, -0, 4, 1.1];
F = sparse([3 3 3 4 4 4 0.0 0.0;
1 2 1 7 9 3 0.0 0.0;
7 8 1 4 8 8 0.0 0.0;
1 1 2 8 5 1 0.0 0.0;
5 6 1 4 2 5 4.6 0.0;
1 6 8 1 1 4 6.8 0.0]);
for i in 1:100
Ik .+= i; # Here is only a simple example of mutation,
F .= F * R .+ sparse(1:6, Nh, Ik);
end
Ah there is a sparse matrix by sparse matrix multiplication happening that was missing in your above examples.
1 Like
Also I am getting the following error:
UndefVarError: I not defined
Can you fix your code above to at least run?
make sure you run the code snippet after posting first. this has happened two times now
Do you mind updating your latest code with a minimal example of simple numbers that still uses sparse matrix multiplicatoin. I would prefer to work with a smaller example, instead of with all the numbers, specific to your use case.
The code still does not run:
ERROR: MethodError: no method matching +(::Vector{Float64}, ::Int64) For element-wise addition, use broadcasting with dot syntax: array .+ scalar
Can you please check again?
Edit: I think you can replace the offending line with Ik .+= 3
Shouldnβt Ik be a vector here, instead of a scalar?
Yes, I am sorry. I have updated it.
You can replace this with just
Nh = [6, 6, 6, 8, 7, 7]; # fixed to reduce the memory use
Thank you for your tip. Any way to get ride of creating the sparse matrix at this line?
F .= F * R .+ sparse(1:6, Nh, Ik);
Instead of creating a dense matrix and then creating a sparse matrix from it, you can directly create a sparse matrix for this line like so:
R = sparse(2:8, 1:7, 1., 8, 8)
I think I have an idea on how to get rid of the for loop even with the matrix multiplication, but it heavily depends on the structure of R
2 Likes
Actually, the for loop is a time-domain loop in my simulation. So, I think, it is kind of mandatory.
What about the structure of R? Is it an off diagonal diagonal matrix?