Purna Chandra Mansingh

Hello, @brandonwillard. Please accept my apologies for forgetting to link my analysis on this; I have now updated the description of this PR with my findings on both `solve` methods....

Can we completely replace the existence `solve` with this Numpy `solve` ?

> @purna135 did you benchmark `cholesky` or `solve`? Yes, I only tested `cholesky`, but since `cholesky` is dependent on `solve`, I thought it would be better to work on `solve`...

Yes, let me test `solve` op.

> @purna135 let's revisit the scipy solve Op and explore the option of add broadcasting there. Ok

Although the benchmark is no longer relevant, I tested both `solve` methods out of curiosity and found that the Numpy implementation of `solve` is **2-3 times** faster than the `vectorized`...

> Aside from that, it seems like those comparisons are only demonstrating the difference between the implementations and/or defaults of `numpy.linalg.solve` and `scipy.linalg.solve` for sparse/broadcasted identity matrices. Does that sound...

I would like to work on this issue

What should be the signature for [Subtensor Op](https://github.com/aesara-devs/aesara/blob/739bd49f996974811592bd3384662d9284426b6d/aesara/tensor/subtensor.py#L682) and [Shape Op](https://github.com/aesara-devs/aesara/blob/739bd49f996974811592bd3384662d9284426b6d/aesara/tensor/shape.py#L43) ?

> If you're talking about constructing symbolic graphs, the signatures are ultimately determined by their `Op.make_node` implementations. Yes, got it now