Ideally, there should be a `random_element`

method but it is not implemented (yet). There seems to be a number of algorithms out there, so if you plan to implement one, please consider contributing it into Sage !

If you do not really mind the density, you can just look at the QR-decomposition of a random real matrix (though this itself is not well defined):

```
sage: G = MatrixSpace(RDF,3)
sage: M = G.random_element()
sage: M
[ -0.1404815835403903 -0.907286820365103 -0.6367014795215493]
[ -0.8102701316983039 -0.6717220161934689 -0.8229332100998983]
[ 0.07523133592237663 -0.12920794641969602 0.8627681516027308]
sage: M.QR()[0]
[ -0.1701173836118859 0.9573741300667837 -0.23344132211739152]
[ -0.9812035951583826 -0.14266646456180654 0.12994531441225488]
[ 0.09110203423256512 0.25115942142054726 0.9636489840135694]
```

Or as a one-liner:

```
sage: MatrixSpace(RDF,3).random_element().QR()[0]
[ -0.7533460901267146 0.23805009753414352 -0.6130267690360484]
[ -0.6227231008312081 0.041418062637191194 0.7813453038052681]
[ 0.21138970689689907 0.9703693601672336 0.117037159341282]
```

If you mind having the density corresponding to the Haar measure, according to https://en.wikipedia.org/wiki/Orthogo... you can generate a real matrix with Gaussian entries such that the diagonal ones are positive:

```
sage: dim = 3
sage: M = Matrix(RDF, dim, dim, [normalvariate(0,1) for _ in range(dim^2)])
sage: for i in range(dim):
M[i,i] = abs(M[i,i])
sage: M
[ 0.0795518237100608 0.07691005207967933 1.0416355962282533]
[ 0.4941624083721146 0.06751228996968504 1.6387527474726973]
[ 1.8797460779164274 0.05562094144225769 0.8965938466096621]
sage: M.QR()[0]
[-0.04089557143145095 -0.8217108088199684 -0.5684354835033859]
[ -0.2540363391790334 -0.5416652758107634 0.8012891284407673]
[ -0.9663297007821968 0.17717244607805865 -0.18659269475761167]
```