StateSpace with all SparseArrays converts to dense on creation

[imciner2]

Trying to create a StateSpace object with all sparse arrays for A, B, C and D leads to the actual StateSpace object having dense matrices instead of sparse ones. If all the matrices given are sparse, then the default statespace type should preserve the sparsity of them. This is extremely important when the system is large and has a lot of zeros in the A matrix.

Reproducing script:

A = spzeros( 2, 2)
A[1,1] = 1
A[2,2] = 1

B = spzeros( 2, 1 )
B[2, 1] = 1

C = spzeros( 1, 2 )
C[1,1] = 1

D = spzeros( 1, 1 )

s = StateSpace( A, B, C, D )
issparse( s.A )

Sample output:

julia> A = spzeros( 2, 2)
2×2 SparseMatrixCSC{Float64,Int64} with 0 stored entries

julia> A[1,1] = 1
       A[2,2] = 1
1

julia> B = spzeros( 2, 1 )
2×1 SparseMatrixCSC{Float64,Int64} with 0 stored entries

julia> B[2, 1] = 1
1

julia> C = spzeros( 1, 2 )
1×2 SparseMatrixCSC{Float64,Int64} with 0 stored entries

julia> C[1,1] = 1
1

julia> D = spzeros( 1, 1 )
1×1 SparseMatrixCSC{Float64,Int64} with 0 stored entries

julia> A
2×2 SparseMatrixCSC{Float64,Int64} with 2 stored entries:
  [1, 1]  =  1.0
  [2, 2]  =  1.0

julia> B
2×1 SparseMatrixCSC{Float64,Int64} with 1 stored entry:
  [2, 1]  =  1.0

julia> C
1×2 SparseMatrixCSC{Float64,Int64} with 1 stored entry:
  [1, 1]  =  1.0

julia> D
1×1 SparseMatrixCSC{Float64,Int64} with 0 stored entries

julia> s = StateSpace( A, B, C, D )
StateSpace{Float64,Array{Float64,2}}
A =
 1.0  0.0
 0.0  1.0
B =
 0.0
 1.0
C =
 1.0  0.0
D =
 0.0

Continuous-time state-space model

julia> s.A
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

julia> issparse( s.A )
false

[olof3]

I had a look at the code and there is indeed a problem in the constructor for the StateSpace type.

You could try using HeteroStateSpace which allows arbitrary matrix types for all the four matrices, although it is not as well supported as StateSpace.

There is some ongoing work on the type system. Hopefully there will be some more progress on this after the summer.

Would be interesting to know more about your use case. Are things too slow when you are using dense matrices and what operations is it that are problematic?

[mfalt]

Did you try the explicit constructor?

StateSpace{Float64,SparseMatrixCSC{Float64,Int64}}(A,B,C,D,0.0)

I think that should work as a work-around.

Yes, hopefully we have time to finish the rewrite of the type system soon. We have been discussing if there are applications warranting having StateSpace{T,<:AbstractMatrix{T}} or if it we should simplify it to be mäMatrix and rely on HeteroStateSpace for the case with sparse matrices.

One problem with we have encountered is that there is not always a good constructor for a general Matrix type SomeMatrix{elType}(data), when we want to promote to a specific matrix with a specific element type.

Edit: it might be that this line have to be changed for my suggestion above to work.

ControlSystems.jl/src/types/StateSpace.jl

Line 55 in d183cb1

return StateSpace{T,Matrix{T}}(MT(to_matrix(T, A)), MT(to_matrix(T, B)), MT(to_matrix(T, C)), MT(D), Float64(Ts))

To
StateSpace{T,MT}(MT(to_matrix(T....

Anyway, this is definetely something that we will make sure works without any complicated work-around once we finish up PR #295

[olof3]

It seems to work if you use the input exactly as suggested by @mfalt .

However, if you just put 0 (an Int) for the sampletime, then it doesn't work due to the problematic line pointed out in the previous post. But, yeah, this should be sorted out properly.

[imciner2]

While I understand the problems with finding the overarching matrix constructor, I don't think splitting the StateSpace type between dense and sparse is a useful thing to do. I think more types of StateSpace should be focused more on how the actual systems differ from each other rather than the storage of those systems in memory. For instance, it would make sense to eventually introduce a DescriptorStateSpace type (so it includes the mass matrix of the system) - and if the normal statespace type is split across matrix types then that one would need to be as well (meaning 2 new types were introduced instead of just 1).

[baggepinnen]

The matrix type is no longer a parameter of the regular state-space type, the following works for sparse matrices

julia> HeteroStateSpace(A,B,C,D)
HeteroStateSpace{Continuous, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}}
A = 
 1.0  0.0
 0.0  1.0
B = 
 0.0
 1.0
C = 
 1.0  0.0
D = 
 0.0

Continuous-time state-space model

julia> HeteroStateSpace(A,B,C,D).A
2×2 SparseMatrixCSC{Float64, Int64} with 2 stored entries:
 1.0   ⋅ 
  ⋅   1.0