Optimizing matrix multiplication algorithm
I have to admit I'm a bit confused, you did the hard parts and have a question on the easy parts? I might be misunderstanding your question.

I am also keen on understanding how I can make it robust such that it doesn't crash when non-rectangular matrices are passed as argument i.e. matrix that has uneven number of elements in different Vds it holds.

You could validate that you have a well formed matrix by
inline bool isValid(const Md&amp; mat)
{
    if (mat.m.size())
    {
        int size = mat.m[0].v.size();
        for (int i = 1; i &lt; mat.m.size(); i++) {
            if (size != mat.m[i].v.size())
            {
                return false;
            }
        }    
    }
    
    return true;
}


and incorporating it into the matrixDot function for validation similar to the shape validation you have right now
    if (m1.m[0].v.size() &amp;&amp; m2.m.size())
        if (m1.m[0].v.size() != m2.m.size())
        {
            LOG(&quot;Shape mismatch: &quot; &lt;&lt; &quot;matrix1 columns: &quot; &lt;&lt; m1.m[0].v.size() &lt;&lt; &quot;, &quot; &lt;&lt; &quot;matrix2 rows: &quot; &lt;&lt; m2.m.size());
            throw std::exception();
        }
        
    if (!isValid(m1))
    {
        LOG(&quot;Invalid matrix :: &quot; &lt;&lt; std::endl);
        std::cout &lt;&lt; m1;
        throw std::exception();
    }
    
    if (!isValid(m2))
    {
        LOG(&quot;Invalid matrix :: &quot; &lt;&lt; std::endl);
        std::cout &lt;&lt; m2;
        throw std::exception();
    }

Any optimizations I can think of are around utilizing std::array instead of std::vector as you are aware of row and column lengths at time of creation.